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^^ Abstract 

A wave-ice interaction model for the marginal ice zone (MIZ) is reported, which involves both the 

I I attenuation of ocean surface waves by sea ice and the concomitant breaking of the ice by waves. It is 

,S^ specifically designed to embed wave-ice interactions in an operational ice/ocean model for the first time. 

Wh We investigate different methods of including the wave forcing, and different criteria for determining if 

/-\ they cause floes to break. We also investigate and discuss the effects of using various attenuation models, 

rrt finding that predicted MIZ widths are quite sensitive to the choice of model. Additional sensitivity tests are 

• performed on: (i) different parameterizations of the floe size distribution (FSD), including the initial FSD 

fi used; (ii) the properties of the wave field; and (iii) the sea ice properties such as concentration, thickness 

• "T^ and breaking strain. Results are relatively insensitive to FSD parameterization but vary noticeably and 

f^.,^ systematically with its initial configuration, as they do with properties (ii-iii). An additional, somewhat 

r^ surprising sensitivity, is the degree of influence of the numerical scheme that performs wave attenuation 

rS I and advection. This is because a naive implementation of spatial and temporal discretizations can cause 

' ' the waves to be over-attenuated, leading to a reduction of the predicted MIZ width by a substantial factor. 

Example simulations intended to represent conditions in the Fram Strait in 2007, which exploit reanalyzed 

h^^ wave and ice model data, are shown to conclude the results section. These compare favorably to estimates 

I of MIZ width using concentrations obtained from remotely-sensed passive microwave images. 

00 

fS| 1 Introduction 

>— ^ In this paper we present a parameterization of wave-ice interactions for an improved representation of marginal 

rv^ ice zone (MIZ) dynamics. The work is part of a project aimed at enhancing forecasting abilities by creating the 

T-H first wave-ice component for use in an operational ice/ocean model. The demand for predictions of wave-ice 

lj_ interactions in forecasting models is becoming greater as accessibility to the ice-covered regions of the ocean 

•^H increases due to the impact of climate change (Stephenson et at, 2011). Indeed, Prinsenherg and Peterson 

/\ (2011) recorded flexural failure induced by swell propagating within multiyear pack ice during the summer of 

3 2009 at very large distances from the ice edge in the Beaufort Sea. While the local sea ice there qualified as 

being heavily decayed by melting [Barber et al., 2009), and thus more fragile, these observations suggest that 

such events could occur more frequently well within the ice pack in a warmer Arctic. The coupling will be 

embedded in a HYbrid Coordinate Ocean Model (HYCOM) , initially in the Fram Strait but intended for future 

implementation in more general ice/ocean models and oceanic general circulation models. The main effects that 

we include in our model are the attenuation of waves as they travel from the open ocean into the MIZ, and the 

subsequent strain-induced fracture of the ice caused by the ice-coupled waves. 

While the notion and importance of integrating wave-ice interactions into an ice/ocean model was broached 
more than two decades ago the third author (VAS), it is only recently that a prototype analysis has been 
published [Dumont et al., 2011). That work provides the platform for this investigation. The model used 
by Dumont et al. (2011) is one-dimensional, i.e. a transect of the ocean, yet integration into an operational 
model clearly requires an extension to two-dimensions, i.e. the ocean surface. A few different themes have been 
identified for investigation before this geometrical restriction is tackled. The first theme is the manner in which 



the wave forcing is included and, in particular, the assumption that there is no interaction between waves with 
differing frequencies. The second is the floe-breaking criterion employed, and the third is the sensitivity of the 
model results to the choice of attenuation model. We also perform a sensitivity study on the effects of different 
parameterizations of the floe size distribution (FSD), of properties of the wave field, and of sea ice properties such 
as concentration, thickness and breaking strain. These parameters can cause significant, although systematic, 
variation in the model results, which highlights the importance of having accurate estimates for the values to 
use, and also the need for more experimental observations to be carried out. 

The general study concludes with some numerical considerations like how the choice of the advection- 
attenuation scheme (AAS), the grid size and the time step used, affect our predictions, before we present some 
more realistic simulations in the Fram Strait using wave and ice model data from 2007. The MIZ widths 
predicted by these simulations compare quite favorably to estimates made from concentrations obtained from 
remotely-sensed passive microwave images (AMSR-E, University of Bremen). 

The energy carried by a traveling wave is distributed over a spectrum of frequencies, and it is consistent to 
interpret the effect of the wave forcing on the ice by integrating over the full spectrum. This approach allows 
for constructive (and destructive) interference between frequencies, and we shall call it the integrated spectrum 
method (ISM). The ISM was used in ice-breaking contexts by Langhorne et al. (2001) and Vaughan and Squire 
(2011), and also by Rattier (1992) who used it to estimate rates of fioe collisions and the resulting (acoustic) 
noise production. 

Langhorne et al. (2001) estimated the lifetimes of ice-sheets in the presence of waves, using the model of 
Fox and Squire (1994) to define a relationship between significant wave and strain amplitudes in the ice. This 
is similar to the approach taken in our study, although we prefer to consider a distribution of strains rather 
than a single representative strain value. Vaughan and Squire (2011), on the other hand, describe a method of 
estimating the width of the MIZ according to the attenuation coefficient of Squire et al. (2009) combined with 
a strain threshold for ice fracture. This work also has similarities to the current study, particularly in regard to 
the use of a parametric probability threshold. However, their model docs not output the FSD, and consequently 
is not detailed enough to include in an ocean model. 

Dumont et al. (2011) did not consider interactions between wave frequencies, based partly on the resulting 
algebraic simplifications and partly on the assertion that different wave groups would separate due to dispersive 
effects (as the group velocity, i.e. the speed that the energy propagates at, depends on frequency). We will refer 
to this approach as the wave group method (WGM). 

The computational loads required for the ISM and the WGM are similar. We find in our results section that 
their model predictions are also quite similar when only breaking resulting from strain failure is allowed for. 
(Stress induced failure was also considered in the WGM setting by Dumont et al. (2011); see below.) Moreover, 
when our waves-in-ice model (WIM) is extended to two-dimensions by considering a full directional spectrum, 
we find there is little difference in complexity between using ISM or WGM. However, because the ISM takes a 
more orthodox theoretical approach than the WGM, and also because it produces a smoother FSD, we advocate 
its use in future simulations. 

Having advected our waves into the ice and allowed for some attenuation to occur as they travel in from the 
ice edge, the predicted distribution of wave energy must be mapped to a fracture probability using floe-breaking 
criteria. Dumont et al. (2011) propose two floe-breaking criteria, one based on stress failure and one based on 
strain failure. The former approximates the horizontal stress in the ice, assuming rigid ice and a monochromatic 
wave surface, and was intended to allow for the effects of cavitation and wetting. During cavitation, the ice 
floe does not follow the wave profllc exactly and potentially causes a strong localized stress on the fioe. The 
criterion predicts greater stress when the waves are longer, which is unphysical in this regime as ice is relatively 
unaffected by long waves because of their low slope/curvature, normally small amplitude, and the low velocities 
they force surface objects to move at. As long waves experience the least attenuation from the presence of ice 
cover, the stress criterion results in an unphysically wide MIZ. Consequently, in due course we will eliminate the 
stress criterion from the parameterization. However, in the future it is worth investigating a different method 
of estimating the likelihood of events like cavitation and wetting occurring, as we find that in general only 
allowing for strain failure causes an under-prediction of MIZ widths. Nonetheless such events are more likely to 
be caused by higher-frequency, shorter waves, so any model including them will need to reflect this. 

One difficulty with using strain failure to determine the probability of floe-breaking is that experimental 
measurements for strain thresholds in sea ice are scarce, whereas there exists an extensive database on its 
flexural strength {Timco and Weeks, 2010). We therefore present a method for converting flexural strength to 
strain thresholds. To do this, we must also estimate associated ice mechanical properties such as the Young's 
modulus and Poisson's ratio, which have been measured less often — especially at the strain rates relevant to 
this paper. The flnal result is a strain threshold parameterized by the brine volume fraction, which is a standard 



parameter in operational sea ice models. 

The primary effect of the ice-cover on a traveling wave is to reduce its intensity exponentially with distance. 
A relationship between this phenomenon and wave scattering by inhomogeneities in the ice cover, especially at 
floe edges, has been established. The first model of wave attenuation in the ice-covered ocean was produced by 
Wadhams (1973). Recent appreciation of the role of multiple scattering and wave coherence in the attenuation 
phenomenon has generated a number of more sophisticated models {Kohout and Meylan, 2008; Squire et al, 
2009; Bennetts and Squire, 2012). The result is that there is now a selection of attenuation models available 
for implementation in one-dimensional wave-ice interaction theory, which can provide an attenuation rate as a 
function of the properties of the incident wave field and ice cover. In this paper we compare the MIZ-width 
predictions produced by three different attenuation models and focus on the FSDs they are generated from. The 
first model is calculated by ensemble averaging using a realistic power law distribution for the FSD. However, 
this process produces artificial resonances that allow a window of high-frequency (and thus high-curvature) 
waves to be transmitted with negligible attenuation through the MIZ. The result is that the MIZ becomes far 
too wide. The second model uses the method of Bennetts and Squire (2012), which calculates the average over 
all possible wave phases throughout the domain. This model produces a more realistic level of attenuation but, 
compared to the most thorough attenuation transect conducted to date (collected by Squire and Moore (1980) 
in the Bering Sea on compliant sea ice floes less than 0.5m thick), it still underestimates the attenuation of 
long waves. Therefore the third model adds some additional damping empirically in order to attenuate these 
long waves in accord with how they are observed to diminish. Both the second and third models predict more 
realistic MIZ widths than the first, but there is still a factor of 2 difference between them. This highlights 
the importance of resolving the problem of how long waves should be attenuated theoretically, and also of 
conducting more experiments to confirm and extend the observations of Squire and Moore (1980) to different 
thickness ranges. 

2 Model components 

The interactions between ocean waves and sea ice that we are interested in occur on small to medium scales, 
but they have a profound effect on the large-scale dynamics and thermodynamics of the sea ice. On a large 
scale the ice cover deforms in response to stresses imposed by winds and currents, and thin ice categories 
(< 300 to 500 mm) are rapidly embedded into thicker ice through convergence, ice formation and ridging. It is 
customary to model the ice cover as a uniform viscous-plastic material, with all complexities subsumed into the 
internal stresses of the associated rheology (e.g. Hibler, 1979; Hunke and Dukowicz, 1997). The viscous-plastic 
description and most historical ice model developments are based on observations of the central pack where the 
ice is compact. It is not representative of the part of the ice cover in contact with the open ocean (the MIZ), in 
which the presence of waves produces smaller fioes and lower surface concentrations. The ice in the MIZ also 
provides less resistance to external forcing and air-ocean heat fluxes and is modulated by the longer endurance 
of thin ice types. Separate rheologies for the MIZ have been proposed [Shen et at, 1986; Feltham, 2005), as well 
as parameterizations for floe size-dependent thcrmodynamical processes {Steele et al., 1989; Steele, 1992). The 
HYCOM model, in which the wave-ice component will be embedded, uses wave-ice interactions to delineate the 
MIZ from the inner ice pack. The collisional rheology of Shen et al. (1987), in which the collisions are implicitly 
driven by wave action, is applied in the MIZ and a standard viscous-plastic rheology is used for the remainder 
of the ice cover. 

The bridge between the two scales is the floe size distribution (FSD), which is parametrized according to 
the observations of Toyota et al. (2006, 2011), as discussed in §2.1. Although the viscous-plastic rheology is 
independent of the floe sizes, the FSD plays a key role at the large scale in determining which rheology to use 
at a particular point. Correspondingly, the MIZ rheology, has an explicit dependence on the FSD. 

At the medium scale, waves are advected from the ocean and are attenuated by the ice. The waves also 
induce stresses and strains in the ice that may or may not cause them to break. Breakage alters the FSD, 
lowering the mean floe length and increasing the amount of attenuation, which is predominantly determined by 
small scale scattering events at the edges of individual floes. An important aspect of the (numerical) breaking 
process is the means by which we decide whether breakage will occur at a particular point in space, and this is 
discussed in §3. It is also strongly influenced by the sea ice properties we use — notably, the Young's modulus, 
Poisson's ratio, and flexural strength of the ice, which generally depend on the brine volume fraction (see §2.2) 
— and potentially the way we deal with the issue of fatigue (see §2.3). The attenuation models we use are 
discussed in §2.4. 



Dc = ( ^^ :7^-^ ) , (1) 



2.1 Floe size distribution 

In this section we review recent progress towards understanding the FSD in the ice-covered oceans, and describe 
how this is used in the wave-ice component of the modeL Details of the probabihty density functions (PDFs) 
considered in the investigation are given in Appendix B. Prior to 2006, numerous researchers (e.g. Weeks et ai, 
1980; Rothrock and Thorndike, 1984; Matsushita, 1985; Holt and Martin, 2001; Toyota and Enomoto, 2002) 
had made observations of floe sizes in Arctic areas. It was found that the FSD generally obeyed a power law 
(Pareto) distribution, having a PDF given by (41) with maximum floe size Di -^ oo. However, many of the 
observations were theoretically problematic as the fitted exponent 7 was often greater than 2, which implies an 
infinite ice area if small floes obey this law {Rothrock and Thorndike, 1984). This can also be seen by letting 
Dq — > in (41), in equation (42b). However, we note that this can also be avoided if the power law is not used 
for very small floes (i.e. by keeping Dq non-zero). By a similar argument, 7 should be greater than 2 for larger 
floes (let Di — > 00 in equation 42b to see this) . Furthermore, if Dq = the average floe diameter is also infinite 
if 7 > 1, and the probability of obtaining a floe size greater than zero (which should be 1) is infinite if 7 > 0. 

To address these problems, Toyota et al. (2006), using data obtained from the southern Sea of Okhotsk, 
investigated the FSD for floes ranging between Im and 1.5 km in diameter. They found that there was a 
decrease in the value of 7 (to about 1.15) when the floe size dropped below about 40 m. The regime shift was 
also observed in Antarctica in the late winter of 2006 and 2007 by Toyota et al. (2011), based on observations in 
the northwestern Weddell Sea and off Wilkes Land (around 64°S, 117°E) with a helicopter-borne digital video 
camera. Concurrent ice thickness measurements were also made, using a helicopter-borne electromagnetic (EM) 
sensor above the Weddell Sea and a video system off Wilkes land. The regime shift was consistent with the 
value 

,48pwatcr5(l-J^^ 

which corresponds to the length below which flexural failure cannot occur {Mellor, 1986). (The parameters Y, 
h, V, Pwater and g are the Young's modulus, thickness and Poisson's ratio for sea ice, the density of the water, 
and the acceleration due to gravity.) Toyota et al. (2011) also developed a fractal model that would produce 
the distribution (41) for smaller floes, based on a fragility parameter 

n = C-^ 7 = 2 + log^H, (2) 

which is the probability that a floe will break into ^^ pieces. Note that this model predicts 7 < 2. A value of 
7 > 2 corresponds to H > 1, which implies that there must be additional breaking mechanisms other than the 
simple breaking scheme used to derive the small floe distribution. Since the regime shift occurs roughly around 
diameter Dc, a possible reason is that the smaller floes are mainly produced by breaking due to flexural waves, 
with the exception of very small bits of rubble and brash ice formed from floe-floe collisions. Larger floes are 
primarily influenced by other processes such as the stresses induced by winds and currents. In addition, large 
floes may also be the result of smaller floes recombining in periods of low temperature, and low wind and wave 
activity. 

A method for using wave- ice interactions to calculate the FSD was proposed by Dumont et al. (2011). They 
calculated the maximum floe size possible in a given wave field, -Dmax, and assumed that in each grid cell the 
floe lengths obeyed a single-regime power law (41), setting Di ~ -Dmax and taking Dq to be a set minimum 
floe length Dmin- The resulting distribution was coupled to the fractal model of Toyota et al. (2011) in their 
floe breaking simulations, but only once the value of -Dmax had been determined. Finally, the expected value 
of D, (-D), was obtained using a novel formula derived as part of the investigation. An alternative means of 
calculating (D) would have been to use equation (42a). In their simulations, they used the values Dmin = 20 m, 
^ = 2, and H = 0.9. 

Here, we choose to modify this approach, by using the observed regime change between different power law 
exponents that occurs around Dc- To do this we must modify the coupling between the wave-ice interactions 
and the FSD. If the wave fleld is large enough that breaking occurs, -Dmax reduces, so we increase the probability 
of having small floes, i.e. we increase the parameter Pq — P{D < Dc) in (43) using the values for 71, 72, Dq and 
-Di below, so that ¥{D < D^^x) — Pmax is a high probability that we set to 0.95. (Pq can be found analytically.) 
Thus -Dmax is no longer an absolute maximum but is now a diameter that most of the fioe diameters are smaller 
than. We then calculate the expected floe diameter {D) from (43), using -Dq = -Dmin = 20 m, Di = D^ and 
7i = 1.15 (which corresponds to fragility H = 0.55), since this was observed in both the Sea of Okhutsk {Toyota 
et al., 2006) and in Antarctica {Toyota et al., 2011). For the larger floes, we take the exponent 72 to be between 
2.1 and 3. We do not rule out the possibility that breaking produces -Dmax < Dc, but once this happens we do 



not allow Dmax to decrease any further. Furthermore, -Dmax can never decrease below -Dmin- If ^max < ^cj 
it is possible that we need Pq = 1 in order that ¥{D < Dmax) — Pmax- If this happens, we must then reduce 
Di in (43) from Dc, and Di is then an absolute maximum. An additional complication that may occur is that 
-Dmax is SO large that there are no small floes and Pq = 0. In that case we would need to increase Di from Dc, 
but instead we assume that Umax > -Dunif, where Z?unif is chosen to be a sufficiently large diameter; 200 m is 
used here. In this case all the floes in the grid cell are taken to be of length Dmax- We shall refer to D^^if as 
the FSD cutoff length. 

Another alternative method of changing the FSD in response to breaking is to make Pq dependent upon the 
other parameters in the distribution (43), e.g. by requiring that ps be continuous at D = Dc (which it seems to 
be in the plots shown by Toyota et al. (2011)), and increase the fragility 11, as Toyota et al. (2011) noted that 
this increased closer to the ice edge. This might make an interesting investigation for the future, but as it is 
slightly more complicated than the previous methods because it cannot be done analytically, we do not attempt 
to do this here. 

It was hoped that by using the observed distribution (43), we would implicitly account for some effects 
other than breaking that are present in the MIZ, such as refreezing. However, it turns out that the MIZ widths 
predicted are very insensitive to changing between FSD (43) and the one used by Dumont et al. (2011). Having 
said this, some sensitivity to the initial value of Dmax used does exist. This is tested in §4.2.4 and §4.3. 

2.2 Ice properties 

Timco and Weeks (2010) provide an excellent up-to-date synopsis of our current knowledge about the engineering 
properties of sea ice. This includes information on flexural strength. Young's modulus and Poisson's ratio, which 
will be used for the present investigation. We note that there are many reported measurements for the flexural 
strength and Young's modulus of first-year ice, but far fewer for multi-year ice. Only a very small data set 
is available for the wave-induced failure strain of sea ice in situ. The data that do exist suggest a value of 
approximately (5 ± 3) x 10~^. In contrast, many flexural strength tests have been made due to its importance 
and use in ice engineering problems. Although we recognize that ice fracture is considerably more complicated, 
we will utilize the experimental findings on flexural strength to parameterize the probability of ice breaking due 
to the strains caused by wave loading. Sea ice is a heterogeneous aggregate of solid ice, brine, air and solid 
salts, and thus flexural strength is not a basic material property, but an indicial strength, as non-uniform stress 
fields are created under flexural loads that require assumptions about material behavior. 

Timco and O'Brien (1994) collate and analyse 939 flexural strength measurements conducted by 14 different 
investigators under a variety of conditions and test types, namely, in situ cantilever tests and simple beam tests 
with 3- or 4-point loading, to show that the flexural strength Uc has a very simple dependence on brine volume 
fraction vi,: 

(Jc = o-Q exp (-5.88\/uh) , (3) 

where ctq — 1.76 MPa. This is plotted in Figure 1(a), and shows a monotonic decrease from ctq as wt, increases. 
Brine volume is often a parameter in ice-ocean models, but it can also be calculated from the formula of 
Frankenstein and Gardner (1967): 
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lib = 10-^-;^^ 0.532 + 49.185 x-f^, (4) 

where Sjcc is the salinity (psu) of the ice and Tjco is its temperature (°C), while sP.^ = Ipsu and 2]"^, = 1°C are 
the unit values of each. 

Flexural strength tests are normally analyzed by means of Euler-Bernoulli beam theory, in which the stress 
normal to the beam cross section, i.e. along the beam axis, cth, is related to the strain, i.e. beam relative 
extension, en, by the one-dimensional Hooke's law cth = Yen, where Y is the Young's modulus. Further 
details are given in Appendix A and by Meirovitcti (2010). In principle, therefore, to convert flexural strength 
into a breaking strain Ec for a beam of sea ice, all we need is Y for sea ice, as £c = o-c/Y. However, an ice floe 
is generally too wide to be modeled as a beam, and is usually modeled as a thin elastic plate, so we must also 
modify the breaking strain by a factor involving Poisson's ratio ly (also see Appendix A): e^ — '^c/{Y{l — i^^))- 

Unfortunately, as Timco and Weeks (2010) point out, the constitutive relation for ice is rate-dependent 
and includes a delayed elastic (anelastic) response (primary, recoverable creep), e"^, viscous strain (secondary 
creep) and strain due to cracking (tertiary creep), as well as instantaneous elasticity, e*. For modest strains and 
timescales, e' and e'^ are of particular importance, as both are impermanent although recovery from delayed 
elasticity is not instantaneous. In the course of a typical flexural strength test and during the recurring cyclic 
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Figure 1: Behavior of the flexural strength (a), and our models for the effective Young's modulus (b) and the 
breaking strain (c) with the brine volume fraction Vh. 



flexure imparted by ocean surface gravity waves, it is expected that the sea ice will experience stress levels 
and rates such that the total recoverable strain e^ « e* + e''. This suggests a variation on the instantaneous 
elastic Young's modulus Y that allows for delayed elasticity to act, that is often called the effective modulus 
or the strain modulus and that we shall denote by Y* . The Young's modulus itself, Y, can only be measured 
dynamically, e.g. ultrasonically with small isolated samples or in situ by measuring the propagation of high 
frequency elastic waves, while Y* < Y arises from assuming an elastic constitutive relation in an experiment 
where some recoverable primary creep occurs. Moreover, as with flexural strength, brine volume affects both Y 
and Y* . Timco and Weeks (2010) report a linear relationship for Y{vb) of the form Y = Yo{l — 3.51vb), where 
Yq « lOGPa is roughly the value for freshwater ice at high loading rates. But, whilst increased brine volume 
leads to a reduction in the effective modulus Y* , the data are too scattered for an empirical relationship for 
Y*(vi,) to be expressed. For "average" brine volumes ranging from 50 to lOOppt (wf, — 0.05 to 0.1, Frankenstein 
and Gardner, 1967), this suggests Y will reduce to between ^ 6-8 GPa. 

As we have noted above, the effect of brine volume on Y* is more difficult to pin down, but we believe the 
same kind of reduction would not be unreasonable. More challenging is determining the effect of anelasticity 
(delayed elasticity) on reducing Y to Y*. The mechanisms that achieve this power-law primary creep with no 
microcracking cause relaxation processes to occur during cyclical loading, so the rate of loading is important. 
Few data can help us here but Figure 4 in Cole (1998) shows model predictions for the effective modulus at 
four loading frequencies that include those associated with surface gravity wave periods, i.e. lO^^-lO" Hz (or 
0.01 - 1 Hz), and, incidentally, the reduction in Y due to total porosity, i.e. air plus brine. The latter effects are 
comparable in magnitude to the reductions in Y given above; the effect of rate is about 0.5 GPa as wave period 
is changed from 1 s to 10 s, and about 1 GPa from 10 s to 100 s. We therefore consider a reduction of 1 GPa is 
to be reasonable in our model, and in summary 



Y* =ro(l-3.51ub)-lGPa. 



(5) 



This is plotted in Figure 1(b). An appropriate choice of a value for the effective Young's modulus is important 
from the wave modeling perspective too, as the higher Y* becomes the more energy is reflected at each floe 
present and the greater the attenuation experienced by the wave train. However, because the same value of Y* 
is used to convert from flexural stress to failure strain, the analysis is self-consistent. 

The final property we will need to consider is Poisson's ratio. Langleben and Pounder (1963) determined 
it to be i^ = 0.295 ± 0.009 from seismic measurements, so in most wave calculations involving ice (e.g. Fox 
and Squire, 1991) it is simply taken to be 0.3. In Figure 1(c), we plot the breaking strain using this value of 
Poisson's ratio, combined with (3) and (5) in 

(6) 



y*(i-i^2)- 



The breaking strain has a minimum value of approximately 5.3 x 10~^ when Uf, = 0.15 {Y* — 3.8 GPa). The 
value is approximately constant for Vb € [0.1, 0.2]. It shows an increase for both higher and lower brine volumes 
— the less porous ice is predictably stronger, while the more porous ice is more compliant so will be able to 
sustain more bending before breaking. If Vb — 0.05, Sc ~ 7.1 x 10~^ {Y* = 7.2 GPa), while if Vb = 0.1, the 
breaking strain drops to Sc ~ 5.5 x 10~^ (Y* — 5.5 GPa). If Poisson's ratio 1^ = 0, these values fall by about 
9%, so that the minimum is 4.8 x 10^^. 



2.3 Fatigue 

The above discussion of the anelastic response of sea ice does not preclude the possibihty that floes can gradually 
fatigue due to repeated bending imposed by passing waves. Fatigue, whether of the high-cycle type associated 
with elastic behavior and growth of microscopic cracks that eventually reach a critical size for fracture, or low 
cycle fatigue where the stress is sufficient for plastic deformation, is characterized by cumulative damage such 
that materials do not recover when rested, i.e. they behave inelastically as opposed to anelastically. Accordingly, 
the effective modulus approach described above, which includes only fully recoverable elastic deformation, cannot 
accommodate fatigue. There is, however, a suggestion {Langhorne et al., 1998) that an endurance limit, i.e. a 
value of stress for which a material will retain its integrity even when subjected to an infinite number of load 
cycles, exists for sea ice. This value, 0.6, was determined on stationary shore fast sea ice in McMurdo Sound, 
Antarctica. Dumont et al. (2011) therefore reduced their flexural strength by a factor of 0.6. We, on the other 
hand, have chosen not to do this because (i) the ice and wave conditions change rapidly in the MIZ so, while 
a stress greater than 0.6i7c can cause failure in principle, it may still occur at a timescale that is well beyond 
that associated with the local dynamics (recall that the endurance limit is for infinite time) , (ii) fatigue strictly 
negates the use of an effective modulus, as permanent irrecoverable damage is gradually done to the sea ice 
either by the nucleation and propagation of cracks or by secondary and tertiary creep, and (iii) the fast ice data 
of Langhorne et al. (1998) show considerable scatter, which is a common feature of fatigue experiments even 
for simple materials. We rest content, therefore, with the expression for Y* defined in equation (5), noting that 
fatigue can easily be added at a later point if results indicate that it plays a role. 

2.4 Attenuation models 

The operational model selects the rate of exponential attenuation from a pre-calculated look-up table. This 
action is performed for each grid cell and at each time step and depends on the prevailing wave and ice conditions. 
The look-up table describes the functional dependence of the attenuation rate on the properties of the ice cover 
and the wave field. The attenuation model runs independently of the operational model. 

In this section three attenuation models will be discussed. They are all based on linear wave scattering 
theory. The ice is represented by a thin-elastic plate and only responds to fluid motion in flexure, with its 
position on the surface of the ocean fixed. Wave energy may propagate freely through homogeneous regions of 
the ice-covered ocean. Attenuation is not caused by the loss of energy from the fluid-ice system but is due to 
energy being scattered by imperfections in the ice cover. Motion is conflned to two-dimensional transects, i.e. 
one horizontal dimension and one depth dimension. Scattering, and hence attenuation, is produced solely by floe 
edges. Extensions to scattering by other features in the ice cover, such as cracks and pressure ridges, are possible 
(see Bennetts and Squire, 2012) but the restriction to floe edges is sufficient for the present investigation. The 
fluid surface is perforated with a large number of floes and an incident wave originates at one spatial extreme. 
Assuming time-harmonic conditions of a given period, the profile of the wave in the transect is calculated by 
solving the full multiple scattering problem, i.e. the infinite sum of reflections and transmissions of the wave 
between floe edges. Exponential decay is then a product of localization theory, which relies on positional disorder 
and requires proper consideration of wave phases. The reliance on disorder implies the use of an averaging 
approach. The attenuation coefficient, a, is hence calculated as the ensemble average of the attenuation rates 
produced in simulations that are randomly selected from prescribed distributions. It is natural to calculate a 
non-dimensional attenuation coefficient, i.e. per floe, for these types of problem, but this is easily mapped onto 
a dimensional attenuation coefficient, i.e. per meter, for use in the operational model. 

The attenuation models now differ only in the distribution of floes that they use, although they may also 
incorporate certain approximations to facilitate the solution process. The first model considered, denoted A, uses 
a seemingly plausible choice for the distributions. The floe length distribution is based on the split-power law 
discussed in §2.1 and given in (43) with D = (20 m, Dc), 7 = (1.15, 2.5) and Pq chosen to make the distribution 
continuous. It was derived from observations in Antarctic locations. Floe separations are generated from an 
exponential distribution P(G > 5) = exp(— g/(G)), with (G) = {D){c^^ — 1) and ice concentration c — 0.9. 
Toyota et al. (2011) did not measure a floe separation distribution so we picked this distribution arbitrarily. 
The attenuation coefficient is calculated as the average of 100 randomly generated simulations. 

The second and third models are based on the recent work of Bennetts and Squire (2012). Rather than 
considering spatial distributions, Bennetts and Squire (2012) considered the wave phases as random variables 
and averaged over all possibilities, assuming they were uniformly distributed. They argued that the model is 
not intended as a true replica of the MIZ but rather a series of impedance changes that resemble the transition 
a wave experiences when passing between open water and an ice-covered region of the ocean. Consequently, 



detailed predictions about the exact distribution of wave phases cannot be rehed upon, and an assumption of 
uniformity is the simplest possible in the absence of a more realistic model. In this setting the attenuation 
coefficient may be calculated analytically rather than relying on a numerical approximation. The expression for 
the attenuation coefficient can be simplified further if the floes are assumed to be long, so that only the reflection 
produced by a single floe edge is required, and the attenuation coefficient is then given by a = — 21og(l — |i?P), 
where R is the scattering coefficient by a floe edge of the specified thickness. Model B uses attenuation calculated 
in this manner. 

In several previous investigations it has been identified that scattering models under-predict attenuation for 
large periods. This may be a result of the increased influence of unmodeled dissipative mechanisms, although 
it is unclear which, if any, is dominant in this regime. One plausible candidate is that some secondary creep is 
occurring at longer periods where the flexural straining rates are slower. While this issue remains unresolved 
it is useful to add an artificial viscosity to the scattering model to dampen the large period waves. Bennetts 
and Squire (2012) showed that their model could be modified in a simple manner to accommodate viscosity in 
the form of an imaginary term in the fluid-ice coupled dispersion relation. Modifying the attenuation model B 
in this way gives a new model, C, and the value of the damping parameter was set by matching results to the 
most reliable MIZ experimental data available at present {Squire and Moore, 1980). 

Another failing of all the models is that they don't predict the maxima in the attenuation coefficients 
(rollover) observed by, e.g. Wadhams et al. (1988), below wave periods of about 5-6 s. It has been suggested that 
this could come from additional local wave generation, possibly caused by nonlinear wave generation ( Wadhams 
et al., 1986). However, we will not address this issue further in this paper, which unfortunately could mean 
that we are underestimating the extra strains produced by additional wave energy at low periods. The three 
models are summarized in Table 1. Note that all three models include the draft of the floes, which alters floe 
scattering in some circumstances {Williams and Squire, 2010) potentially affecting results signiflcantly. Draft is 
actually computationally important as well, because it removes artiflcial resonances and thus makes averaging 
over thicknesses redundant. Consequently, the neglect of draft is considered to be ill-advised. 



Attenuation Model 


A 


B 


C 


Draft 


Yes 


Yes 


Yes 


Uniform Phases 


No 


Yes 


Yes 


Calculation Method 


Numerical 


Analytical 


Analytical 


Damping 


No 


No 


Yes 



Table 1: Attributes of the different attenuation models used in the paper. Draft refers to whether or not draft 
is included in the scattering model; Uniform Phase refers to whether or not the phases of the waves incident at 
each edge are uniformly distributed or whether they are associated with the floe length distribution; Calculation 
Method refers to whether the ensemble averaging is done numerically, or using an analytic formula; and Damping 
refers to whether additional attenuation is added to the attenuation predicted by the multiple scattering model. 



3 Waves-in-ice models (WIMs) 



In our results section we will compare the effects of three WIMs. These involve different methods of including 
the effects of waves in the MIZ and different combinations of breaking criteria. WIMs 1 and 2 are based on the 
approach of Dumont et al. (2011), and treat different frequencies as independent wave groups that separate out 
as they travel due to dispersion. For direct comparison with Dumont et al. (2011), WIM 1 uses both stress and 
strain breaking criteria, although, as foreshadowed in § 1, we do not advocate using a stress criterion because 
of the unrealistic breaking it predicts when the waves are long. Although WIM 1 is similar to the model of 
Dumont et al. (2011), certain changes have been made to facilitate comparison with WIMs 2 and 3. WIM 2 is 
identical to WIM 1, except that it only uses a strain breaking criterion. 

The final model, WIM 3, has the fundamental difference of allowing interactions between different frequencies. 
This requires the consideration of integrals (approximated numerically) over the entire wave spectrum. This 
approach is similar to those used by Rottier (1992) to deal with fioe collisions and accoustic noise production, 
and Langhorne et al. (1998), Kohout and Meylan (2008) and Vaughan and Squire (2011) who also deal with ice 
break-up. WIM 3 only uses a strain breaking criterion. 

Before we can describe the WIMs, however, we first introduce the approximations we use for the stresses 
and strains in the ice floes. 



3.1 Strain approximation 

Let us assume that a sinusoidal wave profile 

Wicc{x, t) = Acos{hccX - Lot) (7) 

with period T — 27r/a; is traveling in the ice, where A = AW{lij) and the factor W{uj) — fcice/^wator converts 
the wave amplitude for open water A into a wave amplitude for ice-covered water. The quantity fcwatcr — w^/g 
is the infinite depth wave number for open water, and kice its analogue for ice-covered water found from (39) in 
Appendix A. 3. We note that kicc can be approximated by the open water wave number fc^ator = lo^ / 9 in sonic 
circumstances, as was done by Dumont et al. (2011), making W ^ 1. However, the approximation is not used 
in the current paper, as it affects the FSDs predicted by the model. Note that it mainly affects the lengths of 
floes rather than the MIZ width itself. 

Now, using the formula (38) from Appendix A. 2, the maximum strain in the ice can be approximated by 

£max = ^W{u)A = E{u)A, (8) 

where E{uj) is the strain per metre of wave amplitude. The spectral density function for the strain in the ice 
will then be S{(jj)E'^{ijj), where 5* is the wave spectrum that produces the open water displacement. 

3.2 Stress/pressure approximation 

Dumont et al. (2011) hypothesize that cavitation between the ocean surface and the underside of an ice floe 
would result in a pressure loading that could cause the floe to break. The model they envision replaces the 
ice floe with an elastic beam of length L that rests horizontally on a sinusoidal wave profile centered in the 
trough between two adjacent crests, so that L = Aice- Thus breaking only occurs at wave troughs. We, however, 
consider a beam of length L = Aicc/2, supported at the mean ice-water interface (i.e. between troughs and 
crests), so that forces on the beam are maximal (and approximately equal) at both crests and troughs. These 
forces are approximated here in the same way as by Dumont et al. (2011), and are taken to be the average of 
the pressures (integrated over the whole beam) due to buoyancy and gravity [see Figure 4 of Dumont et at, 
2011]. This impHes that 

Q{Lu)A=-{p)gXi,,W{uj)A, (9) 
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where Q(w) is the force per unit beam width per metre of wave amplitude, and (p) = (picc + Pwator)/2. 

The maximum possible value of QA that a beam can endure, which we shall denote by Qc, can be approxi- 
mated by comparing the wave profile to a three point beam-loading setup, where the load acts halfway between 
the supports with a force F. The maximum force that the beam can withstand without breaking is denoted 
by Fc. From Eulcr-BernouUi beam theory (see Schwarz et al., 1981, and Appendix A.l), Fc is related to the 
flexural strength CTc by 

where b is the beam width, h is the beam thickness and Qc — Fc/b. At this point wc make one further adaptation 
to the model of Dumont et al. (2011), as ice is better represented by a plate model than a beam model. In 
Appendix A we discuss why it is more difficult to bend a plate than a beam, and thus a larger stress is needed to 
make the strain exceed a certain breaking strain in a plate than in a beam. Therefore (10) should be modified 
by increasing CTc to a* = crc/(l — v)"^ and the maximum force per unit width that a plate can withstand is 

Q.^'^. (11) 

3.3 Incident waves 

Most wave data arc given parametrically in terms of the significant wave height Hs , the peak period Tm and the 
mean wave direction {Ochi, 1998; WMO, 1998). However, in our one-dimensional simulations we will only use 
Hg and T\i, and take the waves to be traveling directly into the ice. The significant wave height here is related 
to the variance of the sea surface displacement (see below). For the power spectrum of the incident wave we 
invoke the Bretschneider spectral density function (SDF) 5b given by 

SB(..TM.ft)^^^^^oxp(-^), (12) 



The incident wave spectrum will be moved into the ice using the Advection-Attenuation Schemes (AASs) 
described in §3.6, so that we will always be able to determine the spectral density function (SDF) S{lo) at any 
point in the ice field. 

3.4 Waves-in-ice models 1 and 2 

Dumont et al. (2011) argue that, as a result of both attenuation and dispersion, wave groups of different 
frequency will separate out and the effects of each group would be able to be considered in turn (beginning 
with the faster, lower frequency groups). Integrating the SDF S{uj) with respect to a; provides the mean square 
value of the sea surface elevation (w^), and also the mean square wave amplitude (A^), since it is commonly 
assumed that {w"^) = ^ (A^) . The SDF is a function of the frequency that is built from a random phase-amplitude 
assumption (Holthuijsen, 2007). However, using the SDF to produce an amplitude spectrum £/{uj) = (A^(u;))^/^ 
as a function of frequency is not straightforward. One method would be to integrate over a frequency interval 
[uj — ^A(jj,u! + I Aw] . Unfortunately, £/[uj) then depends on the arbitrary choice of Aw, which is not ideal for 
the purpose of comparing different waves-in-ice models. In order to rule out this arbitrariness, Dumont et al. 
(2011) used a simple dimensional analysis where they considered [S] = [^x/'^/ui] to obtain 

^(w) ~ ,^^^^2ojS{uj), (13) 

where we have included the correction for ice-coupled waves. Despite its unconventional approach, this formula 
does appear to provide a reasonable and conservative estimate of the highest waves of a group although it tends 
to overestimate the amplitude of short waves; the ones that are the most rapidly attenuated. An alternative 
formulation based on spectral integrals is presented later to resolve this issue and quantify the uncertainty. 
Comparison between both formulations show similarities well within our remaining uncertainties. 

Having obtained an expression for the expected wave amplitude, the stress and strain amplitudes are calcu- 
lated and, if they exceed CTc or Sc respectively, the ice is broken. (As discussed in §2.3, Dumont et al. (2011) also 
multiplied a^ by a fatigue coefficient fi = 0.6, but this approach is not adopted here.) Thus, the ice is broken if 
either of the following two conditions are satisfied: 

^(w) > A^ = -^ = t^'lu/r V (14a) 

If breaking occurs for a particular frequency, the maximum fioe size allowed in a particular grid cell is changed 
to Aico/2. As noted above, it is questionable whether (14a) should be included at all for long waves, as it implies 
that very long waves need only a very small amplitude to break ice. Such waves have very low curvature so 
ice fioes are easily able to bend and follow their profile, which challenges the fundamental assumptions made in 
deriving (14a). Nonetheless, an alternative method of including the effects of cavitation and wetting, or other 
breaking mechanisms, would be worth considering in the future, but such methods should reflect the fact that 
breaking events are more likely to be caused by shorter waves than longer ones. Despite the misgivings about 
the stress failure condition (14a), for comparision with Dumont et al. (2011) it is useful to include it in WIM 1. 
Summarizing, we now reiterate the three WIMs we will be using. WIMs 1 and 2, like Dumont et al. (2011), 
both consider wave groups individually, using (13). In WIMl, the ice is broken if cither (14b) or (14a) are 
satisfied, while in WIM 2, the ice is broken only if (14b) is satisfied. In the next section WIM 3 is introduced, 
which, like WIM 2, only applies a strain criterion for ice breaking, but it allows different frequencies to interact. 
It turns out that, while less orthodox in its use of wave statistics derived from the spectral density function 
than WIM 3, WIM 2 produces similar results. 

3.5 Waves-in-ice model 3 

In this model we allow for the interaction between waves of different frequencies by integrating over the entire 
wave spectrum. The mean square sea surface elevation, or the variance in the position of a water particle at 
the sea surface, (w^) — rnolw] can be obtained from the formula 

/■OO 

m„[w] = / a;"S'(a;)dcj, (15) 
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where S{uj) is the SDF at a general point in space; either in the open ocean or within the ice after having 
undergone some attenuation. (We will also use the second moment of the SDF, ?7i2, later on.) The significant 
wave height is related to jtiq by Hs ~ A^m^. 

Wave heights generally follow a Rayleigh distribution, for which the probability of a wave amplitude exceeding 
a certain value Ac is approximately 

P(yl> A,) = e-^'/<^'>, (16) 

{Longuet-Higgins, 1952, 1980), where (A^) denotes the mean square amplitude. If the wave spectrum has a 
narrow bandwidth and non-linear effects arc negligible (low wave steepness), then (A^) = 2?Tio[';i'], so 

P(yl> Ac)=e-'^'/2moM^ (17) 

Forristall (1978) showed that some data were not fitted well by (17), while in response Longuet-Higgins (1980) 
showed that these data were fitted well by (16), putting the difference down to non-linear effects. 
The mean square ice displacement is approximately 

/•OO 

(wfj =mo[wicc], m„[wicc] = / a;"S'(a;)W^2(cj)dw. (18) 

In a similar way to the water wave spectrum, the probability of the amplitude A exceeding a certain value Ac is 
P(A > Ac) = exp(— ^^/2TOo['U'ice])- In addition, we can also estimate the number of waves we expect in a time 
interval Ai, A'w, as 

A+ /»^_L,,. 1 

(19) 

( WMO, 1998). More precisely, this is the number of times we can expect a particle to cross its point of mean 
displacement in a downward direction. Note that equations (17) and (19) disagree with the results presented by 
Cartwright and Longuet-Higgins (1956). Presumably they are typographical errors, as Longuet-Higgins (1980) 
uses equation (16). Ny/ also defines a representative wave period for the spectrum 5 at a given point of 
Tw — Ai/iVw, and a representative wavelength (in ice) of Aw determined from (39). Inside the ice cover, Tw 
plays a similar role to Tm, the peak period parameter in the Bretschneider spectrum. 
We can also define an analogous quantity for the strain. Its mean square value is 

/•OO 

(£2)=mo[e], TO„[£]= / u''S{u)E^oj)dij, (20) 

so Pe — P(i?w > £c) ~ cxp(— e^/2?7io[£]), where E-^ is the maximum strain induced in a floe by a passing wave. 
The probability that none of the iVw waves will cause i?w to exceed Sc is then (1 — Pg)^™. Thus to determine 
if the ice will break wc will require that 

1 - (1 - P,)^" > P„ (21) 

or, equivalently, that Pg > Pc = 1 — (1 — Pc)"'^/^"'. If it is, we then decrease Umax to Aw/2 (unless Aw/2 < i^min)- 
The probability threshhold ¥c is flexible, but Pc — 0.5 is a reasonable choice since we only have two possible 
outcomes (breaking or not breaking). If we now define a significant strain amplitude Es — 2y/mo\e], (21) can 
be written in terms of Eg as 




Es>Ec = Sc^J-2/\og{¥c), (22) 

which shows more explicitly how the strain criterion depends on Sc, Pc and Ny^. We will return to the sensitivity 
of our results to these parameters in §4.2.4. 

3.6 Advection and attenuation 

We have two AASs that we will use to move the waves from the open ocean into the ice. As we shall see, 
which one we use turns out to have an extremely large effect on our results. The first one, which we shall label 
AAS 1, is similar to the one introduced by Dumont et al. (2011). The difference is that it advects the spectral 
density function S itself, as opposed to the mean square amplitude £/. The scheme proceeds as follows. Let us 
discretize our space, time and frequency variables using Xj ~ jAx (j = 0, 1 . . .), t„ = nAt (n = 0, 1 . . .), and 
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Wfc = Wo + kAuj k = 0,1, . . . ,N^. We choose ujq and Au so that periods between 2.5 s and 25s are included. 
Then let Sk — a^At be the distance a wave moves each time step. The frequency-dependent group velocity 
Ofc should depend on space also but this proved very difficult to implement successfully due to the very large 
contrast between the wave speeds in ice and water for high frequencies. Then the Courant-Friedrichs-Lewy 
(CFL) number, i.e. the fraction of a grid cell this distance corresponds to, is Ck ~ Sk/Ax. Like the group 
velocity, the CFL number is also frequency-dependent. It should be less than or equal to 1 to avoid numerical 
diffusion due to waves crossing more than one grid cell each time step, and also aliasing between waves traveling 
in opposite directions. As a final notational shorthand let S^ ^ = S{u>k,Xi,tn)- Then the AAS proceeds as 
follows: 

^fe". - Skf + Ck f^nli - S-f) , (23a) 



kj ~ '~'k,j '^k yi^kj-l ^kj 

-t/ - i^y (23b) 

5fe",. =^fe",.exp(-d^-isfe). (23c) 

Here, d^ ■ is the dimensional attenuation coefBcient that depends on the average floe length (-D") and ice 
concentration c, and is updated every time step after breaking. It is a well-known phenomenon that using 
a CFL number that is too much less than 1 can lead to significant numerical damping of waves. However, 
an additional complication in a situation where physical attenuation also occurs is that when C^ < 1 this 
attenuation is numerically magnified. In the case where Ck ~ \, the wave arrives in a grid cell, possibly breaks 
the ice in that cell, and then leaves the grid cell. However if Cfc < 1, if some breaking has been done and the 
attenuation (in the entire cell) increased, the wave still has to cross a proportion of the cell before it is advected 
to the adjacent cell. Thus, due to its slightly slower speed, it has significantly less energy than it would have 
had if it had crossed a whole grid cell each time step. 

To correct for this, in our second advection-attenuation scheme, AAS 2, we introduce a frequency-dependent 
time step t^ ~ Axj a^ > At, so each wave moves a distance Ax each time step (the CFL number Cfe = 1). Now, 
define time indices Ik^m and Uk,m so that if lk,m ^ n < Uk,m, tn G >-^fc,m = {ijn — l)Tk,mTk\. In other words, 
n = lk,m represents the first time t„ moves into the time interval ^k,rm and n = w^.m represents the last time 
it does so. Moreover, because we have defined C^ = 1 for all k, moving from one interval J^fc^m to another is 
equivalent to moving from one grid cell to another. Although the spatial intervals, or grid cells, are more easily 
envisioned than the temporal intervals ,-^k,rm they arc unfortunately more difficult to program as this involves 
keeping track of individual wave packets. The crux of AAS 2 is that just before waves of a particular frequency 
change grid cell, the attenuation coefficient for those frequencies in the next grid cell is fixed until the waves 
leave the grid cell beyond, and so only the waves traveling behind them experience the increased attenuation 
due to any breaking that they do. 

Bearing the above in mind as we return to our algorithm, we first calculate the spectral density at t = 
(to — 1)t/c from 






F;»,-i = 5(wfc,x„(TO-l)rfe) 

where a^^^ now represents the attenuation when t ~ (m— l)Tfc, and ^,™^^ has been calculated already. As 
mentioned above, to stop the waves getting trapped in the cell, the next time the attenuation gets updated is 
at i = TOTfc. We now calculate the unattcnuated spectral density at i = TO-Tfc by advecting the waves using the 
Tfe time step, which is designed to be more (perfectly) efficient. This gives us: 

^r^ = pT-\\- (25) 

We can now interpolate between t = (to, — I)Tfe and t ~ rriTk to find the spectral density at t — tn using: 

Sk = ak {tn - (to - l)Tk) , (26a) 

■yi:k = F^T' + ^ (-^^k - Frk-') , (26b) 



].k ^ ^^ \'^3,k ^ j,k 

Sl,^y;:,exp[^aj:^'sk). (26c) 
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Table 2: Default values for independent model parameters. FSD stands for floe size distribution, and 

V{D < -Dmax) is the proportion of floes taken to have diameters below Umax- If -Dmax > -Dunif , then the floe 

lengths in a cell are all -Dmax- 



Quantity 


Symbol 


Value 


Ice thickness 


h 


2m 


Ice concentration 


c 


0.95 


Water density 


Pw 


1025kgm-3 


Ice density 


Pi 


922.5 kg m-3 


Gravitational acceleration 


9 


9.81 ms'2 


Poisson's ratio 


V 


0.3 


Brine volume fraction 


Vb 


0.1 


Exponent for small floe FSD 


71 


1.15 


Exponent for large floe FSD 


72 


2.5 


Minimum floe size in FSD 


-^min 


20 m 


FSD cut-off length 


-Dunif 


200 m 


Spatial resolution 


Ax 


5 km 


Time step 


At 


400 s 


Minimum wave period 


2vr/cj3o 


2.5s 


Maximum wave period 


27r/wo 


23.8 s 


Spectral resolution 


Aw 


7.5 X 10-2 s-i 


Breaking probability threshold 


Pc 


0.5 


¥{D < D„,ax) 


^ max 


0.95 



If we are using WIM 1 or 3, having S*"^ now lets us calculate £/ from (13) and thus determine the new floe size 
distribution. Similarly, with WIM 3 we can approximate the integrals Tnoiwicc], 'm2[wicc\ and ?Tio[^]: which are 
required to determine if breaking occurs in that context. 

4 Results 

We begin our results section by comparing the different attenuation models. As model A is found to give 
unphysical results, we concentrate on the best two attenuation models and run some idealized simulations with 
constant ice concentration and idealized thickness profiles to test the sensitivity of our results to the various 
model parameters. Finally we run some realistic simulations using 2007 wave and ice model data from the 79°N 
transect of the Fram Strait between 3°E (where the waves are specified) and where the transect hits the island 
Norske 0er (off the east coast of Greenland) at longitude 17.7°W. 

Table 3: Default values for dependent model parameters. FSD stands for floe size distribution. 



Quantity 


Symbol 


Dependencies 


Value 


Flexural strength 


CTc 


Vb 


0.27 GPa 


Effective Young's modulus 


Y* 


Vb 


5.5 GPa 


Breaking strain 


£c 


^c.Y*,v 


4.6 X 10-5 


Fragility 


n 


71 


0.55 


SmaU floe cutoff in FSD 


Dc 


Y*,h,u,p^,g 


54.6 m 



Tables 2, 3, and 4 list the values used for the model parameters. In all simulations, unless otherwise specified, 
we take the wave speeds to be spatially invariant, and for all frequencies to have CFL value Ck = 1. In this 
case AAS 1 and A AS 2 are equivalent. 
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Table 4: Model variables that are updated each time step in the model, 
floes that are small. 



¥{D < Dc) is the proportion of 



Quantity 


Symbol 


Initial Value 


Maximum floe length in grid cell 


-^max 


500 m 


Mean floe length in grid cell 


(D) 


500 m 


Maximum floe length in MIZ 


Dmiz 


Om 


Width of MIZ 


Lmiz 


km 


ViD < Dc) 


Po 





Significant strain 


Es 





Expected number of waves 


Nw 






4.1 Attenuation results 

Figure 2 shows the attenuation coefficients produced by the different models introduced in §2.4, computed for 
two different ice thicknesses. Because the C curves include an empirical inelastic contribution, they produce 
the greatest attenuation for large periods, and by necessity this is more pronounced for thinner ice as data are 
absent for larger thicknesses. Curves corresponding to model A, which are generated using random floe lengths 
obeying the distribution (43) of Toyota et al. (2011), are markedly different from the other curves. Due to the 
small values of average floe length {D) (in Figures 2a and b, {D) is respectively 32.5 m and 52.6 m), there is only 
a very small amount of attenuation of long waves, which qualitatively contradicts the observations of Squire and 
Moore (1980), noting that their data set was for thinner (< 0.5m) Bering Sea ice. There is also some additional 
fine structure in the attenuation from model A for lower periods. In particular, there is an interval of periods 
between about 6s and 12s (the interval moves to higher periods as ice thickness increases), where there is much 
less attenuation than the other models. This will have a profound effect on the floe breaking produced by model 
A, as waves from that range of periods can produce very large strains if they remain unattenuated. 

Figure 3 further illustrates the deflciencies of model A. Figures 3(a) and (b) show how a Bretschneider 
spectrum [Hs = Im, Tm = 7 s) would change after being attenuated by either 250 or 500 floes, using models 
A and C. The attenuated spectra produced by model A have much smaller peak periods than those of C (and 
B), and much larger maximum values. Figures 3(c) and (d) show how both the signflcant wave height H^ and 
the significant strain Es change with N , the number of floes that the waves have passed. After only a small 
number of floes it can be seen that Hs and Es for A (red curve) are several orders of magnitude larger than for 
the other attenuation models, which are roughly the same. We can also see that for model A, Es remains well 
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Figure 2: Behavior of the different attenuation models from §2.4 with period for thicknesses (a) 1 m and (b) 
2 m. The models all include draft, and C contains additional (empirical) attenuation for long waves. For ease 
of reference in Figure 3(a), wave periods of 6s, 7s and 8s are indicated as vertical dotted lines. 

above the approximate breaking strain for the range of values of N that are plotted. Es curves for B and C both 
drop quickly below Ec, suggesting that the width of the MIZ, Lmiz, will be similarly small under either of these 
models but will be significantly larger under model A if strain failure is the main breakage mechanism. In fact, 
in our simulations involving model A, we found that our 450-km transect was almost always entirely broken. 
We therefore disregard model A for the remainder of the numerical results, on the basis that the predicted 
attenuation rates are insufficient to replicate what is observed. 
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Figure 3: (a): Examples of a Bretschneider spectral density when the peak period is 7s (green) and Hg = Im. 
(b-d): attenuation of the incident Bretschneider spectrum from (a), when the ice thickness is 2 m. (b): spectral 
density (SD) after traveling past N — 250 floes (solid curves) and N — 500 floes (dashed curves). Two 
attenuation models are used; A (red) and C (black). The ice thickness is 2m. (c-d): Significant wave height 
and strain amplitude after traveling past N floes, using three attenuation models: A (red), B (solid black) and 
C (solid black). In (d), the strain that Eg must exceed to produce breaking, Ec, is plotted as a dashed green 
line. (Here we have used e^ = 5 x 10"'\ Vc = 0.5 and TVw = 25, so Ec = 3.7 x lO"'''.) 



4.2 Idealized floe breaking experiments 

In this section we run some idealized experiments and observe how the final FSD varies with the choice of WIM 
or attenuation model, the properties of the incident wave spectrum, and the values of the concentration and 
thickness. We also check the sensitivity of our results to the breaking strain and FSD parameterization, the 
wave speeds and whether dispersive effects can be neglected, and to numerical settings like the choice of AAS. 
Sensitivity to the time step and grid size is inferred from changing the wave speeds and the breaking strain. 

The domain in all of our experiments is 500-km-long (consistent with our simulations in the Fram Strait in 
§4.3), and is divided into one- hundred 5-kni-wide cells. Cells 10 to 50 (counted from left to right) contain ice of 
constant concentration, but the thickness is given by a variable profile of the form 



hj = 



J < 10. 







(27) 



This represents a thickness ramping up from near zero at the ice edge to a constant value h after about 
a;* = 60 X 10'^ m, and is included to simulate in a simple way the thicknesses at the ice edge obtained from 
the model data from the TOPAZ reanalysis {Sakov et at, submitted) that is used in §4.3. While not fitting 
Fram Strait observations precisely, the thicknesses of sea ice floes in MIZs generally do ramp up from relatively 
modest values near the ice edge, where the ice may be melting or being destroyed by rough seas, to larger 
thicknesses in the interior. Such profiles have been observed both in the Antarctic {Worby et at, 1988) and in 
the Greenland Sea {Wadhams, 1983) from north of Fram Strait down to the Denmark Strait. 

The left hand grid cell (cell 1) contains wave energy fitting a Bretschneider SDF (12) that propagates to the 
right into the ice, breaking it as it goes. The wave field in the left hand cell is kept constant in time, so that 
the final floe length distribution is in equilibrium for this sea state. 

4.2.1 Sensitivity to the choice of attenuation model and the waves-in-ice model 

Figure 4 shows the floe length distribution as we move into the ice cover from the open water, with the different 
rows corresponding to when the ice is broken by WIMs 1, 2 and 3 respectively. The two columns have different 
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Figure 4: Length distributions produced by the different attenuation models and WIMs. The WIMs used are 
WIM f (a, b), WfM2 (c,d) and WfM3 (e, f). The soUd curves correspond to attenuation model B, while the 
dashed curves correspond to attenuation model C. The significant wave height used in the Bretschneider SDF 
(f2) is 3m, while the peak periods arc 7s (a,c,e) and 8s (b,d,f). The thickness parameter h is 2m, and the 
concentration has a constant value of 0.75. 



incident wave spectra; the left column has Tm = 7s and the right has Tm = 8s. Like Dumont et al. (201f), 
we also define two parameters to summarize the floe sizes in the MIZ and the MIZ width. In the context of 
our results, we define the MIZ as the region of broken ice between the open ocean and the unbroken ice, where 
Dmax climbs to its initial setting of 500 m. As Figure 4 shows, this climb is very sharp. The width of the MIZ 
we shall call Lmiz, and we shall label the maximum fioe size in the MIZ Dmiz- For example, in Figure 4(c) 
attenuation model B produces Lmiz ~ 60 km and -Dmiz ~ 98 m. 

Figures 4(a, b) use WIMl, which is the model that is most similar to that used by Dumont et al. (2011). 
Figures 4(c,d) use WIM 2, in which the stress condition has been removed. We observe that this shortens the 
MIZ considerably. Figures 4(e,f) use WIMS, in which the ISM is applied. The results are similar to those 
produced by WIM 2. 

One general aesthetic difference between the results of WIMs 2 and 3, which is also seen in Figure 5, is that 
the -Dmax curves are smoother using WIM 3 since the breaking length is estimated from the entire spectrum. 
With WIM 2 (and WIMl), the same curves are step-like, since we only consider each frequency separately in 
the WGM method. Thus the breaking length there depends more upon the spectral discretization that we use. 
This smoothness of the WIM 3 results will be numerically advantageous when trying to couple with a large-scale 
ice rheology. Wider MIZs could also result from using WIMS, caused by strains due to different frequencies 
interfering with one another, both constructively and destructively, potentially leading to the possibility of more 
breaking for a given wave spectrum. However, it seems that the opposite is actually observed, as both Figure 4 
and Figure 5 show. This is possibly due to the way the wave statistics are obtained from the SDF in the 
WGM. As mentioned earlier, (13) overestimates the mean square amplitude for lower periods — the ones that 
produce more strain and thus do the breaking in WIM 2. Nevertheless, the difference is not significant. There 
are also very small differences between the values of Dmiz, which are slightly longer in WIMS. This is because 
the WGM simply picks the shortest wavelength that will break the ice, whereas the ISM considers the whole 
wave spectrum, and longer wavelengths that might not necessarily be able to break the ice by themselves still 
contribute to the integrals defining the expected value of A\v and thus pull the expected wavelength up. 

When the attenuation model is changed from B to C we can see that using B generally results in wider 
MIZs. This is especially apparent for the higher value of the peak period Tm, as might be expected since the 
added damping in the attenuation is designed to be stronger for longer waves. 

Figure 5 shows the behavior of the MIZ's maximum fioe size, Dmiz, and its width, Lmiz, as we vary the 
peak period Tm continuously while using the different WIMs and attenuation models with some different ice 
thicknesses and significant wave heights. This figure shows that there is a clear difference between attenuation 
models B and C, with B resulting in much wider MIZs for the larger peak periods. This is particularly true 
for WIMs 2 and S, and confirms the importance of including the empirically-based damping for long periods. 
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Figure 5: The behavior of Dmiz and Lmiz with the Bretschneider peak period and significant wave height, Tm 
and Hg, the attenuation model and WIM used. The black, blue and red curves correspond to Hs = f.5m, 
2.5 m and 3.5 m respectively, while solid and dashed curves correspond to the attenuation models B and C 
respectively. WIM f is used in (a, b), WIM 2 in (c, d) and WIM 3 in (e, f). The ice thickness is 2 m and the 
concentration is 0.75. 



However, this also highlights the importance of both explaining theoretically the extra attenuation that is 
observed, and of obtaining more attenuation measurements for large wave periods to confirm the results of 
Squire and Moore (1980). 

When we vary the significant wave height Hg, by looking at the transition from the black curves to blue to 
red, we see that, unsurprisingly, Lmiz increases as we increase Hg. This is true for both attenuation models B 
and C. Similarly, Lmiz also increases with Tm, due to longer waves being attenuated less. I?miz also shows a 
general increase with Tm (for both B and C), but it shows a lot more fine structure than Lmiz — particularly for 
WIMs 1 and 2. The steps in those curves are due to the spectral resolution, as the wavelengths corresponding 
to low periods are relatively far apart. For example, when lo = u}2 {T — 15.82 s), Aico/2 ~ 180 m, but when 
uj — ll)3 {T — 12.36s), Aicc/2 ~ 136m. In addition, in Figure 5(a), we can see in the blue attenuation model A 
curve (WIM 1, Hg — 2.5 m) that -Dmiz unexpectedly decreases to 136 m after having already increased to 180 m. 
This is because the ice is completely broken and the restriction on the distance means that ^/{u!^) has not been 
able to be attenuated enough to prevent it breaking the ice in the last grid cell. If the ice cover was made long 
enough then I?miz would increase monotonically as expected. 

Comparing the overall results of the different WIMs, we see that the MIZ widths are again much larger for 
WIM 1 than the others. For smaller peak periods, Dmiz is also larger when we use WIM 1. Again WIMs 2 and 
3 are extremely similar, so from here on we will just assume they give approximately the same results and only 
compare WIMs 1 and 3. Similarly we will assume that attenuation model C is more reliable than B but we will 
bear in mind that B produces wider MIZs. 

4.2.2 Sensitivity to ice thickness and concentration 

As we increase the ice thickness we expect the MIZ width Lmiz to decrease due to the greater attenuation 
associated with stiffer ice that doesn't bend so readily. This is indeed observed in Figure 6 when we compare 
the dashed curves in 6(b,d) (ft, = Im) with the solid ones {h — 3m). As seen in earlier plots, WIM 1 produces 
a wider MIZ than WIM 3 because of the added effect of breakup induced directly by stress condition (14a) in 
the former case. For this model, Lmiz increases monotonically for both thicknesses with the peak wave period 
until it reaches 450 km, i.e. until all the ice is broken. The effect of significant wave height Hg is predictable — 
higher waves cause more breakup and an increase in Lmiz- The MIZ widths produced by WIM 3 are about half 
the WIM 1 widths. 

Figures 6(a, c) also show that, in general, the maximum MIZ floe length Dmiz depends strongly upon the 
thickness h, approximately doubling as h increases from 1 m to 3 m for both models. As in Figure 5, the WIM 1 
curves exhibit a step-like structure due to the spectral resolution. However, they still generally increase within 
a monotonically increasing envelope. When WIM 3 is used, in Figure 6(c), the curves all behave in a simple 
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Figure 6: How Lmiz a-nd -Dmiz vary with ice thickness; attenuation model C with an ice concentration of 0.75 
is used throughout. Black, blue and red curves respectively correspond to Hg = 1.5 m, 2.5 m and 3.5 m, while 
the dashed and solid curves designate ft, = 1 m and /i = 3 m respectively. In (a, b) WIM 1 is used while in (c, d) 
WIM 3 is used. 



manner, increasing smoothly and monotonically. This is also what was observed in Figure 5. The significant 
wave height makes little difference to Dmiz for either model. 

In general, an increase in concentration increases the attenuation and so causes the MIZ width £miz to 
decrease. This is seen in Figure 7, where the dashed curves correspond to c = 0.5 and the solid ones to c = 0.95. 
All curves increase monotonically with Tm and the wider MIZ that follows from a bigger Hg is predictable. The 
-Dmiz curves in Figure 7(c) do not appear to be affected very much by the concentration, but Lmiz roughly 
doubles when the concentration is dropped to 0.5. 




Figure 7: The effect of a change of concentration on Lmiz and Dmiz for attenuation model C with the ice 
thickness set at 2 m. The black, blue and red curves correspond to Hg = 1.5 m, 2.5 m and 3.5 m respectively, 
while the dashed and solid curves designate c = 0.5 and c — 0.95 respectively. In (a, b) WIMl is used while in 
(c,d) WIM 3 is used. 
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4.2.3 Numerical considerations 

There are a number of numerical approximations we can make during the implementation of this model. For 
example we would like to know the effect of the type of advection-attenuation scheme (AAS) we use, and the 
sensitivity to the time step At and the grid size Aa;. In addition, we will investigate whether it is admissible 
to approximate the wavelengths for waves in ice by the wavelengths in water, and whether it is justifiable to 
ignore the effects of dispersion and set the wave speed constant for all frequencies. 

In results not shown, we found that approximating the wavelengths by values for water tended to increase 
the width of the MIZ only slightly, but reduced the value of Dmiz by as much as half. Although our primary 
interest in the current investigation is in Lmiz , the reduced floe sizes could potentially affect other aspects of a 
coupled ice-ocean model (e.g. thermodynamics), and therefore the approximation is considered ill-advised. 

However, when testing the effects of dispersion we did replace the spatially variable group velocity for waves 
in ice by the open water group velocity. This was because using a spatially variable wave speed to allow for 
different ice thicknesses in different grid cells proved to be numerically problematic, primarily because of the 
large difference between speeds in open water and ice for small wave periods. Specifically, the CFL number 
of these lower-period waves Cj^k was as low as 0.03 while they were traveling in open water, and they were 
thus numerically attenuated almost entirely before they reached the ice. However, since low-period waves are 
attenuated the most strongly, we would not anticipate that approximating the group velocity in this way would 
have a large effect, and so the results obtained will still give us a good idea of how dispersion affects our MIZ 
widths and the floe sizes inside the MIZ. 




Figure 8: The effect of changing the AAS on Lmiz and -Dmiz- WIM3 is used with attenuation model C, the 
concentration is set to 0.75 and Hg = 3 m. The ice thickness is 2 m in (a, b) and 3 m in (c, d). The solid lines 
correspond to AAS 2 and the dashed lines to AAS 1. The time step is set to a constant 260 s and wave speed is 
set to CAx/At for all wave frequencies, where C = 1 (black curves), 0.9 (blue curves), 0.8 (red curves) and 0.7 
(green curves). 



In Figure 8, we begin by looking at the advection scheme. There we set the wave speed constant for all 
frequencies and the time step equal to a constant At = 260 s. Physically speaking, varying the speed should 
thus have no effect on -Dmiz and Lmiz-, and it is comforting that the results produced by AAS 2 show only 
small (numerical) variations for the different CFL numbers used, as illustrated by the close grouping of the 
solid curves. The solid black curve (which agrees with the dashed black curve, i.e. both schemes agree when 
C = 1) is the benchmark to compare the other curves against. However, the results produced by AAS I show 
a dramatic drop in MIZ width as the CFL number C is decreased. Even using C — 0.9 (dashed blue curve) 
consistently reduces Lmiz by a half, and C = 0.7 (dashed green curves) reduces it by three-quarters. These 
results are consistent for both the ice thicknesses shown, h — 2m. (8 a, 8 b) and /i = 3 m (8 c, 8 d). The values of 
Z^Miz are also smaller, but the effects are not as dramatic as they are for Lmiz- As discussed earlier in §3.6, the 
explanation for the difference is that in AAS 1, waves are overly attenuated for C < 1, as the breaking they do 
is effectively moved ahead of them (since Dmax and (D) are uniform over the whole grid cell) , and they must 
then move through this broken ice to 'escape' into the next cell. AAS 2 corrects for this effect by delaying the 
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update of the attenuation coefRcients until the waves move into the next cell. 

Considering Figure 8 now from a physical point of view, and assuming the AAS 2 results are the correct ones, 
we see that the width of the MIZ increases monotonically with the peak period of the incident wave spectrum. 
This is because, although the initial significant wave height is the same, the small wave periods are the most 
attenuated (see Figure 2, for example), so Hg decays more slowly into the MIZ as Tm increases. This generally 
corresponds to larger strains and more breaking further into the MIZ. However, since the dominant wave period 
Tw also increases with both Tm and penetration, and longer waves produce less strain, we would expect that 
the -Lmiz curves would reach a maximum and drop off further if Tm increases to much larger values. 



200 r 




Figure 9: Same as Figure 8, except that the effect of dispersion is accounted for, with waves of frequency w 
traveling with their group velocity in open water, but scaled so that the maximum speed for waves with periods 
between 2.5s and 25s is {CAx/At). 

In Figure 9, we observe the results produced when we introduce the effects of dispersion. In this case the 
results of both AASs are largely insensitive to C, the CFL number of the fastest wave. The wave speeds are the 
group velocities of waves in open water, but scaled so that the fastest wave has speed CAx/At. The global time 
step is again set to 260 s for consistency. The Lmiz results from AAS 1 are very similar to the results shown in 
Figure 8 for the lowest CFL number. In general the AAS 2 results for Tmiz sue only slightly smaller than their 
Figure 8 counterparts, with the exception of two local minima at about Tm — 7.5 s and 8.5 s when /i = 2 m (9 b) 
and at Tm ~ 7.5s and 8.2s when h — 3ni {9d). The minima are much less pronounced for the larger thickness. 
Although they drop suspiciously close to the AAS 1 results, they are very robust with respect to changes in the 
CFL number and appear to be genuine physical phenomena. 

To explain the presence of the minima in AAS 2 let us reconsider Figure 8. There, the leading wave- front 
contains waves of all frequencies traveling at the same speeds. The leading waves have only traveled through 
unbroken ice and are therefore relatively unattenuated. Consequently, a particular grid cell experiences a few 
time steps of large waves passing it before the more attenuated waves arrive and Hg starts to decrease. In 
addition, it is during these few time steps when Hs is at a maximum that breaking will happen if the waves are 
large enough and, if breaking docs happen, it causes additional attenuation which decreases Hg even further. 
When dispersion is introduced (as it is in Figure 9), each grid cell will also experience this general pattern of Hg 
building to a maximum and then dropping again as time progresses. Also, as before, the maximum Hg and the 
dominant wave period Tw will increase with Tm. However, the leading wave front will only consist of the longer, 
faster waves and the maximum significant wave height will be less than if dispersion was neglected. Moreover, 
the shorter waves produce more strain and they, being slower, will have also had to travel through broken ice, 
so less breaking will occur further in. This explains why Tmiz is smaller with dispersion than without. Also, 
because Hg is smaller in each cell the probability of breaking occurring is much closer to the critical probability. 
Therefore even though a slightly higher incident peak period produces higher significant wave heights further in, 
because Tw is also higher there may also be slightly lower strains and so the breaking probability may drop to 
just below Pc. And, of course, increasing Tm further increases Hg and so the breaking probability may increase 
again. 
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Figures 9(a,c), and a couple of the AAS 1 curves in Figures 8(a,c), also show a step-like structure in Dmiz, 
similar to what was observed in Figures 5-7 for WIMsl and 2. In that case, the steps were due to the 
spectral resolution, but this is not such an issue here since we are integrating the wave spectrum. The steps in 
Figures 9(a,c) and Figures 8(a,c) can better be attributed to the same mechanism that we proposed to explain 
the minima in the Lmiz curves in Figures 9(b,d). 
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Figure 10: (a-b): The effect of the FSD and the value of the endurance breaking strain e^ on Lmiz- In (a) 
/i = 2 m, while in (b) /i = 3 m. In both plots, when the FSD described in §2.1 is used, Lmiz is plotted as a dashed 
blue curve when e^ = 10^^, a dashed black curve when Sc — 3.25 x 10^^, a red curve when Sc — 5.5 x 10^^, a 
solid blue curve when Sc — 7.25 x 10"^, and a solid black curve when Sc = 10~^. When the FSD is changed 
to the one used by Dumont et al. (2011), Lmiz is plotted for the same values of Ec using the markers V, V, 
A, A and A respectively, (c-d): The effect of Dinit, the initial value of £>max, and the shape of the thickness 
profile on Lmiz- In (c) /i = 2m, while in (d) /i = 3m. Blue, red and black curves correspond to -Dinit = 250m, 
500 m and 750 m respectively. Solid curves correspond to thickness profiles of the form (27), while dashed curves 
correspond to a uniform thickness being used. In all four plots, WIM 3 is used with attenuation model C and 
Ha = 3 m. 



4.2.4 Other sensitivities 

In Figures 10(a, b) we investigate the sensitivity of the MIZ width to the FSD parameters, and the parameters 
involved in the breaking criterion (22) for WIM 3: the breaking strain £c — <Jc/{Y*{1 — v^)), the number of 
waves per time step A'w and the probability threshold Pc. Since Ec is linear in ec^ we expect that it will produce 
the most variability of the latter three. For example, if A'w is taken from the normal distribution A/'(25, 5) and 
Pc from A/'(0.5, 0.15), where Af{m, s) has mean m and standard deviation s, then varying these two parameters 
only produces a standard deviation in Ec of about five or six percent. Also, we note that varying TVw acts as 
a proxy for varying the time step. Therefore, if we let Sc experience a 10-15% range of variability we should 
be able to infer the effects of changing At and Pc. However, although the value of Poisson's ratio i' appears 
to be reasonably consistent for sea ice {Langleben and Pounder, 1963), the value we use for Y* and hence £c 
contains a large degree of uncertainty. As a result, we have let Sc vary over an order of magnitude to obtain a 
more conservative estimate of its effects. 

The curves corresponding to different values of Sc in Figure 10, show that Lmiz is highly sensitive to this 
parameter. This is particularly apparent in Figure 10(a) when h — 2u\. The middle curves in both figures relate 
to Ec = 5.5 X 10^^, the value resulting from using Wf, := 0.1 in (3-6). Clearly reducing £c makes more difference 
than increasing it, but there is still considerable variation for the higher values. Nevertheless, 7.25 x lO^'^ is 
probably close to the upper limit for e^ while 4 x 10^^ is close to the lower limit, so our MIZ widths have an 
uncertainty of about 25% in them due to uncertainties in Eg. Likewise, the variability caused by Ai and Pc is 
about half this. 

Comparing the two FSDs used in Figures 10(a, b), i.e. the curves and their corresponding markers, we can 
immediately see that only negligible differences in the large-scale properties are produced by changing the FSD 
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Figure 11: Model data for our one-dimensional simulations in the Fram Strait in 2007, between the south-east 
coast of Norske 0er (latitude 79°N, longitude 17. 7° W), which corresponds to a: = 438km on our one-dimensional 
grid, and latitude 79°N, longitude 3°E, which corresponds to x = km. The wave field is specified at a; = km 
and is obtained from the WAM ERA-Interim reanalysis. The significant wave heights and peak periods are 
plotted in (a). The waves are then advected west through ice with concentrations and thicknesses taken 
from a TOPAZ reanalysis. They are interpolated onto a regular grid with longitudinal resolution of 0.125° 
(Ax = 2.65 km), and are plotted in (b) and (c). For comparison, the ice edge and edge of the MIZ estimated 
from AMSR-E concentrations are also plotted in (b) as solid and dashed green lines respectively. 



from the one we introduced in §2.1 to the one used by Dumont et al. (2011). Of course, the differences would 
be visible on a smaller scale. It is worth noting that we have treated all floes as experiencing the same stress 
field. In actual fact, the strains in a floe are dependent on its length and shape and accounting for this could 
increase sensitivity to the FSD. However, determining where and how a two-dimensional floe would break in 
such a study is not trivial and will not be attempted here. 

Moving on to Figures 10(c,d) we observe that the initial FSD has more influence. When the thickness 
profile from (27) is used (solid curves), we see that varying the initial value of Umax? which we shall refer to as 
-Dinit, between 250 m and 750 m causes a variation in Lmiz of about 10-25%. There is more variation (25-45%) 
when we use a constant thickness profile, and we note that the MIZ widths are more narrow than their variable 
thickness counterparts. This is because, for the variable profile, waves travel with little attenuation through 
the thin ice near the edge, and the ice is therefore broken. Further in, the thicker ice causes rapid attenuation, 
thus acting like a barrier, and preventing further breakage. In §4.3, when the variation in the ice thickness is 
increased, we shall find that this reduces the sensitivity to -Dinit to around 10% for most of the year. 

4.3 Realistic experiments in the Pram Strait 

In this section we show the results of some simulations using WIM3 with realistic wave forcings, ice concen- 
trations and ice thicknesses for the Fram Strait during 2007. Figure 11(a) shows the wave forcing, which was 
obtained from the ERA-Interim reanalysis (WAM model), while Figures ll(b, c) show ice concentrations and 
thicknesses obtained from a TOPAZ reanalysis {Sakov et al, submitted) in which concentration data derived 
from AMSR-E (University of Bremen) has been assimilated. On average, the modeled ice edge is 45 km west of 
the ice edge observed by AMSR-E, which is plotted as a solid line in 11(b). This is well within the uncertain- 
ties and resolution of the model (TOPAZ has a resolution of about 13 km) and the resolution of the AMSR-E 
analysis (the transect from 15°W to 5°E was divided into bins with widths of about 21.2 km, i.e. 1 degree in 
longitude, and analyzed for ice concentration (Kloster and Sandven, 2011)). The internal concentrations also 
compare well. 

Also plotted (dashed line) in Figure 11(b) is an estimate for the inner edge of the MIZ, determined from the 
same AMSR-E concentrations using the criterion that c < 0.9 corresponds to the MIZ. While this is a different 
criterion from the floe size criterion we use in this paper, it should still give us a rough idea of how well our 
WIMs are performing. 

The mean ice thickness is roughly 0.8 m, creeping up towards 2 m in the summer, which is thicker due to 
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Figure 12: Results of one-dimensional simulations in the Fram Strait in 2007. (a): the lines separating broken 
and unbroken ice, as determined by WIM3 using attenuation model C, are plotted for fi^ = \ (black line), 
/3/i = 1.75 (red line) and fi^ — 2.5 (blue line), where /3/i is a factor used to increase the ice thicknesses from 
Figure 11(c) which are unrealistically low. The breaking strain e^. is 5.5 x 10~^, determined from (6) using 
Vh = 0.1, and the initial value of Mnax in each grid cell is Dinit — 500 m. (b): same as (a) but with e^ reduced 
by 25% to 4.125 x 10^^. (c): same as (a) but using attenuation model B. (d): same as the red curve in (a), but 
using I?init = 250 m (black curve), 500m (red curve) and 750 m (blue curve). For comparison, the edge of the 
MIZ estimated from AMSR-E concentrations is plotted as a dashed green line in all plots. 



greater movement south of multi-year ice from the Arctic Ocean at that time. These ice thicknesses are probably 
too low, so we have also run simulations in which the ice thicknesses are multiplied by a factor /3/j of either 1.75 
or 2.5 in order to get closer to observations (e.g. Widell et ai, 2003). 

Figure 12 shows the results of numerical experiments using the model outlined in this work. Results for 
the expected breaking are calculated using either the TOPAZ thicknesses or the increased thicknesses, the 
TOPAZ concentrations and the ERA-Interim waves. The values of Ec used in each are 5.5 x 10^^ in 12(a, cd) 
and 4.125 x 10~^ in 12(b), and the attenuation model used in 12(a, bd) is either C or B in 12(c). Each night 
in 12(a-c), the floe sizes are re-initialized to being uniformly Z?init=500m long. An extension to including a 
memory in each cell of -Dmax and a gradual refreezing was rejected on the basis that we are unable to include 
the more important effects of ice movement due to winds and current in this one-dimensional experiment. 

In Figure 12(d), the value of I?init is varied to test the model's sensitivity to this parameter. With the 
inclusion of the more realistic variations in ice thickness, this is usually less than 10%. 

In all plots, except for Figure 12(c) when the attenuation model B is used, the edge of the broken ice 
compares well with the AMSR-E-determined MIZ edge in the non-summer months especially, recalling that the 
ice edge from TOPAZ is on average about 45 km too far to the west. In the summer, less breaking occurs due 
to smaller waves and thicker ice. It is also more likely that, due to the ice being more dilute in the summer, 
a concentration criterion for the MIZ will disagree with a floe size criterion. The results presented here show 
that our simulations have a very low sensitivity to the breaking strain, and generally do not vary much with 
thickness, which is one of the least well-known properties. The variations that do occur are systematic in that 
increasing the thickness makes the MIZ narrower, as observed in §4.2.2. The variations with respect to Ec are 
consistent with our previous findings, with smaller values resulting in a wider MIZ. As with the sensitivity tests 
done on Dinit, there is less variation than was observed in the idealized tests in §4.2.4. 

However, Figure 12(c) shows that using attenuation model B produces much more variation with thickness, 
with far too much breaking produced for the lower thicknesses. This re-emphasizes the importance of obtaining 
a better understanding of the attenuation process, and in particular the lower-than-observed attenuation of 
long waves that arises due to inelastic processes. It also highlights the urgent need for more measurements of 
attenuation, and also of more reliable thickness data. 
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5 Conclusions 

There are many observations that suggest a primary role for ocean waves in shaping the morphofogy of ice 
fields, especially in the vicinity of the ice edge where waves are particularly fierce, but also throughout the MIZ 
where they habitually limit the size of the constituent ice floes by fracturing those floes that are too large to 
exist as the waves permeate further into the ice pack. Attenuation, arising due to scattering and supplementary 
inelastic processes such as turbulence, bending hysteresis and interfloe collisions and rafting, also occurs causing 
a gradual reduction of the wave energy envelope with distance from the ice edge that, crEteris paribus, results 
in a gradual increase in floe size with penetration. The FSD is therefore continuously modified by pervasive 
incident ocean wave trains that, according to their period, may either travel long fetches from distant storms or 
else be more locally generated. They are then preferentially filtered by the sea ice in a manner that favors the 
survival of longer wavelengths. Notwithstanding these outcomes, waves can also redistribute ice floes spatially 
by herding them together, potentially creating zones of higher or lower concentration, and also polynyas, and 
they can assist complementary thermodynamic processes to melt the ice cover. 

Given these several influential factors relating to the composition of MIZs, in particular, it is perhaps 
surprising that wave-ice interactions are not yet included in ice/ocean models and oceanic general circulation 
models. While this has been discussed in the past — indeed VAS proposed it in the early eighties — the 
complexity of doing it has proved insuperable until now. The current paper, which provides a road map of 
how to do it, gives some flrst predictions of how floe size and MIZ width are manipulated by waves in a one- 
dimensional spectral setting, and a realistic example of the adapted HYbrid Coordinate Ocean Model ice/ocean 
code in action. It continues the progress made by Dumont et al. (2011) towards embedding a ubiquitous and 
dominant physical process in forecasting ice/ocean models that are used for operational as well as research 
purposes. It sets the scene and provides the machinery to deal with the next stage of development, which 
is to incorporate two-dimensional interactions arising from a directional sea comprising energy at a comb of 
different frequencies distributed angularly and provided cither by observations or a wave model such as WAM 
or WAVEWATCH III®. 

We summarize here the most substantive observations and conclusions wc have reached during the reported 
research programme. 

1. Three models describing how waves break up the sea ice were tested; WIMs 1 and 2 use individual wave 
groups, while WIMS allows wave frequencies to interact. WIMs 2 and 3 break ice when the strain Ec is 
exceeded, while WIMl also applies a stress condition for direct comparison with Dumont et al. (2011). 
WIMs 2 and 3 produce very similar results, while WIM 1 generally produces MIZs that are too wide. 
This is because long waves in WIM 1, which are relatively unattenuated, are (unphysically) more likely to 
produce breaking than short waves. 

2. Various options exist to describe the way ocean waves attenuate as they propagate through open ice 
fields, differing in the type of fioe size distribution used, the averaging process employed and any imposed 
approximation to facilitate solution. Three models are used in this paper; they are labeled A, B and 
C. Models B and C are anchored in the recent work of Bennetts and Squire (2012) who recognized that 
accurate information about wave phases could not be provided by two-dimensional scattering models 
(one horizontal and one vertical dimension), and so phases should simply be assumed to be uniformly 
distributed. Although it was derived from a supposedly realistic FSD, model A did not do this with 
the result that artificial resonances were created that produced the over-transmission of high-curvature 
waves and thus too much breaking. Model C also included an empirical parameterization of wave energy 
loss from inelastic phenomena other than scattering. This additional damping made a large difference to 
our results, giving MIZ widths that were roughly half those produced by model B. This highlights the 
importance of resolving the problem of how long waves are attenuated theoretically, and also of conducting 
more experiments to confirm the observations of Squire and Moore (1980) and to extend them to different 
thickness ranges. 

3. Various sensitivity tests were done using the different WIMs. It was found that a larger significant wave 
height gives a wider MIZ, which is common sense, and that the higher the ice concentration the narrower 
is the MIZ and the smaller is the maximum fioe size. Increasing the ice thickness also decreased Lmiz, 
but increased -Dmiz • MIZ width increased monotonically with the peak period of the incident spectrum. 

4. It was found that different parametcrizations of the FSD had only a negligible effect on MIZ width and 
maximum floe size. The initial FSD had a noticeable influence on the idealized results of §4.2 but less 
influence in the realistic simulations of §4.3. 
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5. The results in §4.2 were also found to depend strongly on the failure strain e^ of sea ice. While very few 
measurements of £c are available, a large data set has been compiled for flexural strength that is consistent 
with Sc at the strain rates concerned. Nevertheless it still contains large uncertainties and should be 
investigated further experimentally. Again, however, the realistic simulations showed less sensitivity to 
this parameter. 

6. Two advection and attenuation schemes (AASs) are used as described in §3.6, AAS 2 being more physically 
realistic. The less sophisticated AAS 1 was found to produce too great a dependence of MIZ width on 
wave speed, while AAS 2 performed well. 

7. The realistic simulations in Fram Strait, 2007 compare reassuringly well with contemporaneous AMSR-E 
(University of Bremen) data. They were conducted using the WAM ERA-Intcrim reanalysis to provide 
the wave forcing and a TOPAZ reanalysis for ice concentrations and thicknesses. 

A Elastic beam and plate results 

Ice flexural strength tests are typically analysed using Eulcr-Bcrnoulli beam theory (EBBT) {Meirovitch, 2010). 
Suppose the neutral axis of an unloaded Euler-BernouUi beam runs along the x (xi) axis, and suppose the z (x^) 
axis points upwards, in the opposite direction to the loading direction. Also let the beam have width b, thickness 
h and length L. If z S [—h/2, h/2] and x € [— L/2, L/2], the EBBT approximates the strain en as en = zd^w, 
where w{x) is the vertical deflection of the beam, and relates the stress an to Sn by the constitutive relation 
(Til = Yen, where Y is the Young's (elastic) modulus. The other property of the beam that we need to apply 
EBBT is the second moment of area, given by 

b/2 i.h/2 ,73 

/ z^dzdy^^. (28) 

-b/2 J -h/2 ^^ 

The equation of motion for an Euler-BcrnouUi beam equation is 

dl {YIdlw{x,t)) ^ q{x,t) - piccbhdfw{x,t), (29) 

where q is the load per unit length and p\cc is the ice density. (For a static formulation, such as during a flexural 
strength measurement, q{x,t) — q{x) and w{x,t) — w{x), so any time derivatives of w disappear.) If Y and / 
are constant, then (29) becomes 

^dtw{x, t) - ^^ - pMM^, t), (30) 

where the forcing q/b is now a pressure. 

In Appendix A.l we solve this equation in one of the three most common testing situations to find the 
stresses and strains in the beam and thus also the breaking strain and fiexural strength. (We do this for the 
three-point-loaded beam test). The expression for the fiexural strength (34) derived in this way is found to be 
the same as given by Schwarz et al. (1981), and the associated breaking strains allow us to use the wealth of 
flexural strength measurements to inform our strain breaking criteria. This is not a new result, but we think it 
is helpful to have a derivation written down in one easily accessible place as Schwarz et al. (1981) neither derive 
their expressions nor give any references which explain them. 

After this we present formulae for the stresses and strain in an Euler-BcrnouUi thin elastic plate, which is 
a more accurate model for an ice floe in the MIZ than the beam model. The essential difference between the 
two models is that for a beam, the width b is small enough that Poisson's effect is not important — the top 
of the beam is under compression, so the width there increases slightly, while the bottom of the beam is in 
tension so the width shrinks — thus on average the width stays roughly the same. For a plate, however, the 
width is of the same order as the length, so this cancellation is not present to the same extent, and strains in 
the y direction must be accounted for. This effectively makes a plate stiffer than a beam, as to bend it in one 
direction you must bend it in the perpendicular direction as well. This can also be seen by comparing (30) to 
the Eulcr-Bcrnoulli thin elastic plate equation 

Ace V''w(x, y, t) = p{x, y, t) - picehd^wix, y, t), (31) 
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where p{x,y,t) is an applied pressure, V — (V^)^, V^ = 9^ + dy, and Ace = (Fft,^/12)/(1 — v'^) is the flexural 
rigidity {Fung, 1965). The quantity v < 1 therein is Poisson's ratio and the 1 — v^ term increases the stiffness 
of the plate by roughly 10% (for v k, 0.3) from what it was in the beam equation [Yh^ /12). 

An important consequence of this is that the same stresses produce smaller strains in a plate than in a beam, 
making a plate harder to break. In other words, even though the stress may exceed the flexural strength for a 
beam, the strain may not exceed the breaking strain (which we assume is the same for an ice beam and an ice 
plate) and the plate will not break. 

A.l Simple beam loaded at three points 

Here a beam of length L is simply supported at its ends x — ±i/2 and loaded in the centre {x — 0) with a 
downward force F until it breaks. (When this happens F — Fc.) The beam displacement w satisfies 

YIdiw{x) ^ -FS{x), (32a) 

w{±L/2) = 0, (32b) 

YIdlw{±L/2) = 0. (32c) 

Note that the delta function has units of m^^ since its integral with respect to a; is 1. Now, since the general 
solution to d'^f{x) = d{x) is f{x) = -\x\ +Ax + B, where A and B are arbitrary constants, we can immediately 
write d'^w as 

dlw{x)^^^{L~2\x\). (33) 

(An equivalent way of formulating the equation d'^f{x) ~ S(x) is to require that d'^f(x) = for |a:;| > 0, and 
that lim \dxf{x + e) — dxf{x — e)] == 1.) We can integrate (33) further and apply (32b) to find w, but we do 

not need to as we are only interested in the strain. The strain is given by 

Fz 
e{x)^zdlw{x)^^^^{L-2\x\), 

which has a maximum absolute value of Emax = {FhL)/(8YI) when a; = and z — ±h/2. Multiplying this by 
Y also gives the maximum stress, so the breaking strain and flexural strength are 

/i 2 F,hL F,hL 3F,L 

Sr = —OZwiO) — , C7r = = -. (34) 

A. 2 Stresses and strains in an elastic plate 

From Fung (1965), the stresses in an elastic plate with neutral plane at 2; = are given by 

y^ /£l £12 0\ 

a{x,y,z,t) = ^M^i2 £2 w(x,y,i), (35) 

1 - '^ V 0; 

where /^i = V^ — (1 — 1^)9^, C12 = {l — iy)dxdy, and /^2 = ^^^(l^'^)^^- Fung (1965) then uses the constitutive 
relation 

1 I 3 

/ 92 0^0^ \ 

e{x,y,z,t)^~z\dxdy dl \w{x,y,t). (36) 

V ^vV 

Consequently, the stresses and strains have their maximum magnitudes at the plate surfaces z — ±/i/2. 
Also, if w{x, y, t) does not vary in the y direction (such as when a plane wave is traveling through the plate in 
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to derive the strains as 



a direction parallel to the x axis) the maximum strain and stress (assuming v < 1/2 so that ;^/(l — z^) < 1) are 

£max = max{|£ii|} = -max{|9^w|}, (37a) 

o-max = max{|crii|} = "'"^ , (37b) 

taking the maximum in time over one wave period, and over the x interval we are concerned with. Thus for 
£max to exceed Ec = cTc/Y and the plate to break, we need CTmax to exceed a* — crc/(l ~ t^)'^- Consequently, 
we effectively have a more stringent stress condition for a plate than for a beam (where we would only require 
Cmax > o'c for breaking to occur). 

Finally, we note that for a sinusoidal traveling wave of the form (7), with wavelength Ajcc = 27r/fcico, the 
maximum value the strain will ever take is 

£max = T^^hcc ~ -^TT T^" • (38) 

A. 3 Dispersion relation for a floating elastic plate 

Here we present the dispersion relation for a floating elastic plate floating on water of infinite depth. This 
dispersion relation appears in almost any paper dealing with wave scattering by sea ice, taking the limit as 
water depth goes to infinity if necessary (e.g. Fox and Squire, 1991; Squire et ai, 1995). It assumes the water is 
inviscid and incompressible, that the flow is irrotational, and that the displacement at the ice-water interface is 
small enough that the dynamic and kinematic boundary conditions can be linearized. By coupling the resulting 
surface pressure with the floating elastic plate through the p term in (31), this implies that a traveling wave 
with a surface displacement of the form (7) is only possible if fcicc satisfies the equation 

(Acofcfcc + Pwatcr5 " Picchuj'^)hcc = Pwatcr^^^- (39) 

B Probability density functions 

In the current paper we have used two main probability distributions — the Rayleigh distribution and the power- 
law (Pareto) distribution. Kohout and Meylan (2008) use a Rayleigh distribution to generate the random floe 
lengths and spacings for their ensemble averaging, but our main use for it will be in the area of wave statistics. 
The split (two-regime) Pareto distribution was observed by Toyota et al. (2011) to describe floe diameters in 
the MIZ around Antarctica, and Duniont et al. (2011) used this result to model the floe length distribution once 
they applied their floe breaking criteria to deduce the maximum possible floe length for a given wave field. 

B.l Rayleigh distribution 

The Rayleigh distribution has the PDF 

Pk(A.)4*^^K^)' '"'■•^-"' (40) 

[o, for X < 0, 

where s = (D) \/2/n. It will be useful later to be aware that the absolute values of the extrcma of a continuous 
normally-distributed random variable X with mean zero and variance s^ obey p^{A, s) {Cartwright and Longuet- 
Higgins, 1956; Longuet-Higgins, 1980). For our purposes, A will be wave amplitude and the variance s^ will be 
the integral of the spectral density function S{lj). 

B.2 Power law distribution 

The PDF for the truncated power law distribution is 

p^{D;^,Do,D,)^{D^^ ffi?G[i?o,A], (4^) 

iiDi{Do,D,), 
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where D is the floe length and /3p is chosen so that 

r-Di 



/ ' ppiD)dD = fip{D--^ - D-") = 1. 

J Do 



(The ordinary power law distribution has Di = oo and 7 > 0.) This PDF gives the expected values of floe 
length and area (assuming that floes are approximately circular) of 



{D)= I ' ppiD)DdD 
Jd„ 

= :^(^r-^n, (42a) 

(Area) = 7 / pp{D)D^dD 

4 J Do 

{Dl-'-Dl'-'), (42b) 



7r/3p7 ^r,2-7_ n2-7^ 



4(7 - 2) 
respectively. The split power-law distribution has the PDF 

PSP{D- 7, D,Po) - npAD; li,D„,Di) 

+ {l-Vo)pp{D;-f2,Di,^). (43) 

Apart from the exponents in 7 = (71,72) and the diameters in D = (Do,-Di), there is an additional free 
parameter in this PDF, Pq = V{D < Di). 

Acknowledgments 

This work is part of the project 'Waves-in-Ice Forecasting for Arctic operators' (WIFAR), funded by the Nor- 
wegian Foundation of Research and Total E&P. We wish to thank Aleksey Marchenko, Gareth Hegarty and 
David Leslie for helpful discussions, and Francois Counillon for assistance with the TOPAZ data and related 
discussions. 

References 

Barber, D. G., R. Galley, M. G. Asplin, R. De Abreu, K.-A. Warner, M. Pucko, M. Gupta, S. J. Prinsenberg, 
and S. Julien (2009), Perennial pack ice in the southern Beaufort Sea was not as it appeared in the summer 
of 2009, Geophys. Res. Lett., 5^(L24501), doi:10.1029/2009GL041434. 

Bennetts, L. G., and V. A. Squire (2012), On the calculation of an attenuation coefficient for transects of 
ice-covered ocean, Proc. R. Soc. Land. A, 468(2137), 136-162. 

Cartwright, D. E., and M. S. Longuet-Higgins (1956), The statistical distribution of the maxima of a random 
function, Proc. R. Soc. Lond. A, ^57(1209), 212-232. 

Cole, D. M. (1998), Modeling the cyclic loading response of sea ice. Int. J. Solids Struct., 55(31 — 32), 4067 — 
4075. 

Dumont, D., A. L. Kohout, and L. Bertino (2011), A wave-based model for the marginal ice zone including a 
floe breaking parameterization, J. Geophys. Res., 116{GQAQQl), 12pp, doi:10.1029/2010JC006682. 

Feltham, D. L. (2005), Granular flow in the marginal ice zone, Phil. Trans. R. Soc. Lond. A., 363, 1677 — 1700. 

Forristall, G. Z. (1978), On the statistical distribution of wave heights in a storm, J. Geophys. Res., 83{C5), 
2353-2358. 

Fox, C., and V. A. Squire (1991), Strain in shore fast ice due to incoming ocean waves and swell, J. Geophys. 
Res., 96{C3), 4531-4547. 



28 



Fox, C, and V. A. Squire (1994), On the oblique reflexion and transmission of ocean waves at shore fast sea 
ice, Phil. Trans. R. Soc. Land. A., 347, 185-218. 

Frankenstein, G. E., and R. Gardner (1967), Equations for determining the brine volume of sea ice from -0.5 to 
22.9°C, J. Glaciol, 6, 943—944. 

Fung, Y. (1965), Foundations of Solid Mechanics, Prentice-Hall Inc, Englewood Cliffs, New Jersey. 

Hibler, W. D. (1979), A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815-846. 

Holt, B., and S. Martin (2001), The effect of a storm on the 1992 summer sea ice cover of the Beaufort, Chukchi, 
and East Siberian Seas, J. Geophys. Res., 106{C1), 1017-1032. 

Holthuijsen, L. H. (2007), Waves in Oceanic and Coastal Waters, Cambridge University Press. 

Hunke, E. C, and J. K. Dukowicz (1997), An elastic-viscous-plastic model for sea ice dynamics, J. Phys. 
Oceanogr., 27, 1849-1867. 

Klostcr, K., and S. Sandvcn (2011), Ice motion and ice area flux in the Fram Strait at 79N using ASAR and 
passive microwave for Feb. 2004 - Jul. 2010, Tech. Rep. 322, Nanscn Environmental and Remote Sensing 
Center. 

Kohout, A. L., and M. H. Mcylan (2008), An clastic plate model for wave attenuation and ice floe breaking in 
the marginal ice zone, J. Geophys. Res., ii5(C09016), doi:10.1029/2007JC004,434. 

Langhorne, P. J., V. A. Squire, C. Fox, and T. G. Haskell (1998), Break-up of sea ice by ocean waves, Ann. 
Glaciol., 27, 438-442. 

Langhorne, P. J., V. A. Squire, C. Fox, and T. G. Haskell (2001), Lifetime estimation for a land-fast ice sheet 
subjected to ocean swell, Ann. Glaciol., 33, 333-338. 

Langleben, M. P., and E. R. Pounder (1963), Elastic parameters of sea ice, in Ice and Snow, edited by W. D. 
Kingery, pp. 69-78, MIT Press, USA. 

Longuet-Higgins, M. S. (1952), On the statistical distribution of the heights of sea waves., J. Mar. Res., 11, 
245-266. 

Longuet-Higgins, M. S. (1980), On the distribution of the heights of sea waves: Some effects of nonlincarity and 
flnite band width, J. Geophys. Res., 55(C3), 1519-1523. 

Matsushita, M. (1985), Fractal viewpoint of fracture and accretion, J. Phys. Soc. Jpn., 54(5), 857-860. 

Meirovitch, L. (2010), Fundamentals of Vibration, McGraw-Hill. 

Mellor, M. (1986), The mechanics of sea ice, in The Geophysics of Sea Ice, edited by N. Untersteiner, pp. 
165-182. 

Ochi, M. K. (1998), Ocean Waves. A Stochastic Approach., Cambridge University Press. 

Prinsenberg, S. J., and I. K. Peterson (2011), Observing regional-scale pack-ice decay processes with helicopter- 
borne sensors and moored upward- looking sonars, Ann. Glaciol., 52(57), 35-42. 

Rothrock, D. A., and A. S. Thorndike (1984), Measuring the sea ice floe size distribution, J. Geophys. Res., 
89 (CA), 6477-6486. 

Rottier, P. J. (1992), Floe pair interaction event rates in the marginal ice zone, J. Geophys. Res., 97{C6), 
9391-9400. 

Sakov, P., F. Counillon, L. Bertino, K. A. Lisaster, P. R. Oke, and A. Korablcv (submitted), T0PAZ4: an 
ocean-sea ice data assimilation system for the North Atlantic and Arctic, J. Geophys. Res. 

Schwarz, J., R. Frederking, V. Gavrillo, I. V. Petrov, K.-I. Hirayama, M. Mellors, P. Tryde, and K. Vaudrey 
(1981), Standardized testing methods for measuring mechanical properties of ice. Gold Reg. Sci. Technol., 4, 
245-253. 



29 



Shen, H. H., W. D. Hibler, and M. Lepparanta (1986), On applying granular flow theory to a deforming broken 
ice field, Acta Mech., 63, 143-160. 

Shen, H. H., W. D. Hibler, and M. Lepparanta (1987), The role of floe coUisons in sea ice rheology, J. Geophys. 
Res., 92 {C7), 7085-96. 

Squire, V., and S. C. Moore (1980), Direct measurement of the attentuation of ocean waves by pack ice. Nature, 
283(5745), 365-368. 

Squire, V. A., J. P. Dugan, P. Wadhams, P. J. Rotticr, and A. J. Liu (1995), Of ocean waves and sea ice, Annu. 
Rev. Fluid Mech., 27, 115-168. 

Squire, V. A., G. L. Vaughan, and L. G. Bennetts (2009), Ocean surface wave evolvement in the Arctic Basin, 
Geophys. Res. Lett, 5'6(L22502), doi: 10.1029/2009GL040,676. 

Steele, M. (1992), Sea ice melting and floe geometry in a simple ice-ocean model, J. Geophys. Res., 97{C\\), 
17,729-17,738. 

Steele, M., J. H. Morison, and N. Untersteiner (1989), The partition of air-ice-ocean momentum exchange as a 
function of ice concentration, floe size, and draft, J. Geophys. Res., 94{C9), 12,739-12,750. 

Stephenson, S. R., L. C. Smith, and J. A. Agnew (2011), Divergent long-term trajectories of human access to 
the Arctic, Nature Glim. Ghange, 1, 156-160, doi:10.1038/NCLIMATE1120. 

Timco, G. M. W., and S. O'Brien (1994), Flexural strength equation for sea ice, Gold Reg. Set. Technol., 22, 
285-298. 

Timco, G. M. W., and W. F. Weeks (2010), A review of the engineering properties of sea ice. Gold Reg. Set. 
Technol, 60, 107-129. 

Toyota, T., and H. Enomoto (2002), Analysis of the sea ice floes in the Sea of Okhotsk using ADEOS/AVNIR 
images, in Proc. 16th Int. Symposium on Ice, International Association for Hydro -Environment Engineering 
and Research, pp. 211-217, Duncdin, New Zealand. 

Toyota, T., S. Takatsuji, and M. Nakayama (2006), Characteristics of sea ice floe size distribution in the seasonal 
ice zone, Geophys. Res. Lett., 55(L02616), doi:10.1029/2005GL024556. 

Toyota, T., C. Haas, and T. Tamura (2011), Size distribution and shape properties of relatively small sea-ice 
floes in the Antarctic marginal ice zone in late winter, Deep-Sea Res. II, 55(9-10), 1182-1193. 

Vaughan, G. L., and V. A. Squire (2011), Wave induced fracture probabilities for Arctic sea-ice. Gold Reg. Sci. 
Technol, 67(1-2), 31-36. 

Wadhams, P. (1973), Attenuation of swell by sea ice, J. Geophys. Res., 78(18), 3552-3563. 

Wadhams, P. (1983), Sea ice thickness distribution in Fram Strait, Nature, 305(5930), 108-111. 

Wadhams, P., V. A. Squire, J. A. Ewing, and R. W. Pascal (1986), The effect of the marginal ice zone on the 
directional wave spectrum of the ocean, J. Phys. Oceanogr., 16, 358-376. 

Wadhams, P., V. A. Squire, D. J. Goodman, A. M. Cowan, and S. C. Moore (1988), The attenuation rates of 
ocean waves in the marginal ice zone, J. Geophys. Res., 93{CQ), 6799 - 6818. 

Weeks, W. F., W. B. Tucker, M. Frank, and S. Fungcharoen (1980), Characteristics of surface roughness and 
floe geometry of sea ice over the continental shelves of the Beaufort and Chukchi Seas, in Sea ice Processes 
and Models, edited by R. S. Pritchard, pp. 300-312, University of Washington Press, Seatle. 

Widell, K., S. Osterhus, and T. Gammelsr0d (2003), Sea ice velocity in the Fram Strait monitored by moored 
instruments, Geophys. Res. Lett., 50(L018119), doi:10.1029/2003GL018119. 

Williams, T. D., and V. A. Squire (2010), On the estimation of sea- ice thickness using wave observations, 
Dynam. Atmos. Oceans, -^5(2-3), 215-233. 

Worby, A. P., C. A. Geiger, M. J. Paget, M. L. Van Woert, S. F. Ackley and T. L. DcLiberty (2008), Thickness 
distribution of Antarctic sea ice, J. Geophys. Res., 1 1 3 {CC05S92) , 14pp, doil0.1029/2007JC004254. 

World Meteorological Organization (1998), Guide to Wave Forecasting and Analysis, WMO No. 702, 2 ed. 

30 



